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I-    INTR0D0CTI3N 

A.       GENERAL    PRESENTATION    OP    THE   PHENOMENA 

Great  progress  has  been  made  in  the  static  and  dynamic 
analysis  of  complex  structures  through  the  continued  devel- 
opment of  discrete  element  methods  of  structural  analysis. 
Tremendous  improvements  in  computing  power  have  made 
possible  the  study    of  fully    nonlinear    problems. 

The  analysis  of  the  response  of  a  structure  submerged  in 
a  fluid,  is  severely  complicated  by  the  intrusion  of  signif- 
icant fluid^structure  interaction  effects.  Recently,  the 
development  of  a  variety  of  surface  interaction  approxima- 
tions has  provided  the  means  for  a  more  efficient  analysis 
of  the  coupling  between  the  structure  and  the  surrounding 
fluid. 

Computer  codes  for  structural  analysis  are  well- 
developed  so  that  the  fluid-structure  interaction  is,  for 
the  most  part,  handled  by  adding  new  capabilities  to 
existing   structural   analysis   programs. 

It  is  a  well  known  fact  that  the  primary  threat  to  a 
submerged  structure  is  the  shock  wave  that  results  from  an 
underwater  explosion.  However,  the  complexity  of  the  phys- 
ical phenomena  involved,  along  with  the  difficulties  encoun- 
tered in  obtaining  experimental  results  have  been  serious 
drawbacks  for  the  analysis  of  these  types  of  problems.  But 
there  is  a  definite  need  for  investigations  of  large  defor- 
mations and  buckling  problems  in  a  structure  submitted  to  an 
underwater    explosion. 
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B.       OBJECTIVES 

This  study  deals  with  the  nonlinear  response  of  a 
submerged  cylindrical  shell  to  a  shock  wave.  The  existing 
finite  element  code  EPSA  (Elasto-Plastic  Shell  Analysis) 
[Eef.  1]  which  includes  nonlinear  affects  and  a  surface 
interaction  approximation  was  selected  for  the  study.  In 
order  to  alleviate  the  tedious  interpretation  of  results  at 
points  throughout  the  shell,  PAIRAN-G,  a  color  graphics 
system,  was  used.  PATRAN-3  allows  for  a  global  representa- 
tion of  a  guantity  distribution  over  the  structure  rather 
than  the  discrete  representation  given  by  a  computer  output. 
The  objectives  of  this  study  were  to  merge  the  finite 
element  code  EPSA  with  the  color  graphics  system  PATRAN-G, 
and  to  conduct  an  analysis  of  the  response  of  various  types 
of  cylindrical  shells  to  a  spherical  shock  wave  generated  by 
an    underwater   explosion. 
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II.  THE  EPS A  COMPOTES  PROGRAM 

A.  PRESENTATIOM 

EPSA  (Elasto-Plastic  Shell  Analysis)  [Ref.  1]  is  a 
computer  program  developed  by  Waidlinger  Associates  and 
funded  by  DNA/N AVSEA/ONR  for  the  purpose  of  the  analysis  of 
submerged  stiffened  shells  under  dynamic  loadings.  It  incor- 
porates a  number  of  specific  features  which  are  geared  for  a 
more    efficient    analysis. 

In   particular,    EPSA    allows: 

-The    analysis   of   shells   in    an   acoustic   medium,      subjected   to 

both    low   and  high   frequency    shock   loadings. 

-An    efficient   modelling  of    the   elasto-plastic   behavior 

-The      inclusion    of      large      displacsment      effects    to      analyze 

dynamic  buckling   situations    and   post-buckling   behavior. 

-The    modelling    of  stiffensrs  and   internal   structures. 

-The    fluid-structure   interaction    effect 

The  following  sections  describe  the  equations  of  motion 
for  a  submerged  structural  shell  in  an  infinite  fluid 
sujected  to  a  pressure  loading  and  investigate  the  modelling 
of  the  surrounding  fluid.  The  last  sections  are  devoted  to 
the  finite  element  discretization  as  well  as  to  specific 
features  concerning    EPSA. 

B.  EQOATIOHS    OF    MOTION    FOR    THE   SHELL 

Considering  a  thin  shell  of  thiokness  h  ,  volume  7r  area 
R  submerged  in  an  infinite  fluid  (figure  2.  1)  .  The  shell 
stress  resultants   are  defined    from   the   stress   tensor  by   : 
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h/2 

Jin    /-> 


h/2 

Mij  =  J    <*Lj 


(  dt         and  Mij  =    /     ai-»  t    it 

•h/2'  4/2 

The    distribution      of    strains    (a      ,e  ^e      )       is   assumed 

linear,    the   curvatures  at   mid-height    are    (^it^j^l^    • 


to    be 


«„(Mu) 


N,,  {Hu) 


M|l(Wu) 


HnlNit) 


Figure    2.1        Shell   Stress   State, 


Applying  the  principle  of   virtual    work   gives: 


/. 


■/. 


(s}{6e}dR  =       {p}  (Su}dR 


-     /p  MIS* 


}dR 


(2.1) 


Where 

{u}  =  (u1,u2rW  )T  is  the  displacement   vactor   at   each    point 
(s|  =  (Nn,N22  #Ni2#Hn   ,M22#Mi2  )T       is      the   stress      resultant 
vector 

{e}  =  (e11#e22  ,2ei^kil    ^22»2k12)T       is      the   strain      resultant 
vector 
p  is    the  mass   per   unit  area    for   the    shall. 

Tha    symbol    5     will    refer   10    a   virtual   quantity,      and   the   dot 
or   star    denotes   a   differentiation    with   respect   to    time. 
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The  first  ten  of  equation  (2.1)  represents  the  virtual  work 
of  internal  forces  ,  the  second  represents  the  virtual  work 
of  external  forces  (i.e.  pressure,  point  loading,  etc)  ,  and 
the  third  one  represents  the  contribution  of  inertia  forces 
in  the  virtual  work.  Thus,  this  last  term  expresses  the 
effects   of   dynamic    phenomena    on   the    structure. 

C.       FLUID    HODELING 

In  the  case  of  a  submerged  structure,  the  external 
forces  are  the  pressures  applied  at  the  fluid-structure 
interface. 

ks  the  shock  wave  hits  the  stucture  it  gets  reflected, 
thus  inducing  a  pressure  term  pr  .  In  addition,  the  motions 
of  the  shell  will  also  generate  a  radiated  wave,  with  a 
pressure  contribution   pr^ 

Therefore,  the  pressure  at  the  fluid  structure  interface  is 
the    sum   of   the    incident,    reflected   and   radiated   pressures: 


P       =    Pi    ♦    Pr    *    Pra  (2-2> 

Where 

p     =    Total   dynamic    pressure. 

Pj_  =  Pressure  associated  with  incident  free  field  pressure 
wave. 

pr  =  Reflected  pressure  due  to  the  interaction  of  the  inci- 
dent wave  with  the  structure,  the  structure  being  fixed  and 
rigid. 

pra=      Radiated    pressure   due    to      the   normal   movements     of  the 
surface   of  the    structure     in      the   fluid 
ps=    pr*   prais  called  the   scattered  pressure. 

The   methods    for      getting      the    scattered   pressure      ps  will 
now      be   investigated.  Assuming      a      spherical   wave      in     an 
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acoustic  medium      with  a     sound   speed      c    ,        the   pressure     is 
determined    by   the   well  known   wave   aquation: 


c*y   p      »      Ifg  (2.3) 


and   the      proper    conditions      at   the      boundaries   of      the    fluid 
domain   . 

One  alternative  for  solving  the  previous  problem  would 
be  simply  to  use  a  finite  element  discretization  of  part  of 
the  fluid  domain  ,  imposing  a  raiiation  condition  at  the 
boundary  [Ref.  2].  However,  this  would  require  a  very 
significant  fraction  of  the  computing  effort  that  could  not 
be   devoted   to  the   structural   modelling. 

Therefore,      a   more    efficient   way      of    computing  the   scattered 
pressure  is  required. 

The  Doubly  Asymptotic  Approximation  (DAA)  imparts  upon 
the  structural  model  a  surface  loading  composed  of  incident 
and    scattered   waves. 

In  the  high  freguency  limit,  it  can  be  shown  that  the 
scattered  nodal  force  vector  (Fs}  is  related  to  the  wave 
particle  velocities  normal  to  the  structured  surface  by 
[Ref.    3]    : 

{Fs}    =  [A]{US}  (2.4) 

Where  [ A  ]  is  the  diagonal  matrix  of  nodal  areas  (areas  asso- 
ciated with  each  node)  and  (Us)  is  the  vector  of  nodal  scat- 
tered normal  velocities.  Therefore,  in  this  high  frequency 
case  the  shock  wave  is  simply  a  plane  wave  and  equation 
(2.4)  states  that  the  pressure  is  proportional  to  the  fluid 
velocity. 
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In  the  low  frequency  limit  tha  fluid  structure  interac- 
tion   is   governed   by    the   relation: 

{Fs}    =       Cav]{U*}  (2.5) 

Where    {U *  }=  ^-{O-}    is     the    nodal   normal      acceleration   vector 
s  at     s 

and   [  BL  ]  is   the    added     mass    matrix  computed   in   an    incompres- 
sible  fluid. 

Thus,  in  the  low  frequency  case ,  the  loading  is  due  to 
the  motion  of  the  rigid  structure  in  an  incompressible 
fluid,    a   problem  typically    found    in   hydrodynamics. 

When  a  broader  range  of  frequencies  is  considered,  one 
can  combine  the  two  previous  equations  wi*h  the  differential 
equation   governing    the  scattered    pressure   [Ref.    3],    giving: 


UJ*     {Fs}    ♦     pct^Ji    {Fs}     ■    PciVs)  (2-6) 

Where      {Fs*}    =  JL  {Fs  } 

Defining  the  vector    cf  nodal   scattered   pressures    {ps}    by: 

{Ps}     ■      U1"1     (Fs) 

we   get: 

C»V]{PS}     *    pC[A]{ps}     =    pC[Mv]{Os}  (2.7) 


which  is  the  set  of  equations  that  gives  the  scattered  pres- 
sure at  each  node  of  the  wetted  surface  of  the  shell. 
Equation  (2.7)  gives  exact  results  in  both  the  high  and  low 
frequency  limits.  Thus,  DAA  allows  the  modelling  of  the 
fluid-structure  interaction  via  a  coupled  set  of  differen- 
tial equations  at  the  wetted  nodes  of  the  structure  instead 
of  modelling  the  whole  fluid  with  a  finite  element  grid. 
The   use   of      DAA    will   free   some    memory   space      in   the   computer 
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for  a  better  modelling  of  structural  behavior  while  at  the 
same  time  giving  some  reasonably  accurate  results  for  the 
loading  of  the  structure- 
Therefore,  the  equations  for  the  study  of  underwater 
shocks  will  consist  of  a  coupled  set  of  structural  equations 
that  come  from  the  Principle  of  Virtual  Work  and  of  fluid 
equations  at  the   wetted   nodes   that  :oiae    from    DAA    equations. 

D.       FINITE    ELESENT    PROCEDURE 

1 •      Discret  izat ion 

The  principle  of  virtual  work  is  rewritten  for  the 
structure  introduced   in   section   B  [Bef.    4] 

I   {s}T  (5€}dR   -  I    {p}T    {Su}dR    ♦  /   p {ufcfiu}  dR      =    0  (2.3) 

The  surface  of  the  region  is  covered  by  a  quadrilat- 
eral mesh  of  N  elements,  each  having  an  area  k^  .  The 
previous  integral  can  then  be  replaced  by  a  summation  of 
integrals  over    Aj_  .  : 

{s}     (<5e}lR    =        )      {s}.    {Se}.  dR  (2.9) 


I 


i  =  l    Ai 


The    integrations   over   A  .   ar e     then   performed    by   dividing  the 
element   into  four  regions   Ai       (figure   2.2)     We   have   then: 

f  {s}i   C<5e}.  dR     =         £     {s}.T  {5ek}i  k±  (2.10) 

J  .  k=l 


A.  k=l 


and    therefore   equation    (2.9)    becomes    : 

jtsf  {5e}dR      =    £    £     (s)i   (5ak}i    ^  <2-11> 

^R  i  =  l    k=l 

Using   the   same    procedure,    it   can   also    be   derived    : 
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f{p}T{6u}     dR=    J]  J]     (PJT   C6^   ■$ 

Jr  i=l  k=l 


=     {P}      {«q} 


(2.12) 


/■  N        4 

/p{u}T{5u}    dR=     £    £{u>i    *5Q*L  4a    CM](q}C5q)  (2.13) 

"*  i  =  l  k=l 

Whsre  [M]  is  the  nass  matrix  ,  {P}  is  the  vector  of  exter- 
nally applied  forces,  and  {g}  is  the  vector  of  nodal  normal 
displacements   for  the  structure. 


mi 


m 


■>:■: 


Figure    2.2        Srid   Points   in    EPSA. 

By  definition,  finite  elament  discretization  can 
express  the  displacement  {u}  at  any  point  of  an  slement  as  a 
function  of  the  displacements  at  the  corner  points  of  the 
element,    defined   by    {q}    . 


<U}       =    [H]     {q} 


(2.  14) 


Where  [H]   is  a    6xt2    matrix    of   interpolation    functions. 

Combining  the   derivatives  of    {u}    will   give   the  strain   vector 
{e}  .    In   matrix    form,    equation    (2.14)     gives    : 
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{«}    =      [B]    {q}  (2. 15) 

Where  [B]  is  a  6x12  matrix  functioQ  of  the  geometry  of  the 
element  as  well  as  a  function  of  the  previous  deformations 
in   the   case   of    large   displacements. 

Using  the  previous  result,  equation  (2.9)  is  rewritten  in 
the    following   way   : 


t 


{s}T{6e}    dR*   V   V   {s}^CBk]i{5q}iici»       {F}T  {Sq}  (2.16) 

i=lk=l 


Where    {F}    is  the   internal   force  vector. 

Combining  the  previous  equations  in  equation  (2.8)  ,  the 
principle  of  virtual    work   becomes: 

CM]     {q}  =        V  (    &}±     -       {F}.     >  (2.17) 

i=l 

Therefore,  the  principle  of  virtual  work  has  be  =  n  tranformed 
into  a  system  of  ordinary  differential  equations  which  are 
much    more   amenable    to   numerical   treatment. 

In  EPSA,  each  arbitrarily  shaped  quadrilateral 
element  is  defined  by  four  corner  nodes,  each  having  three 
translational  and  no  rotational  degree  of  freedom.  In  order 
to  represent  the  behavior  in  bending  eight  nodes  not  contig- 
uous with  the  element  are  also  ussd  (figure  2. 3) .  Every 
element  accesses  twelve  nodes  and  has  twenty  degrees  of 
freedom:  three  translational  degrees  at  the  corner  nodes 
and  one  corresponding  to  the  displacement  normal  to  the 
surface  for  each  of  the  eight  exterior  nodes.  The  bending 
behavior  (second  derivative  term)  is  expressed  in  terms  of 
the  nodal  displacements  via  a  finite  difference  technique 
[Raf.    4]. 
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Then  the  nodal  displacement  vector  of  element  i  is 
simply: 

(q)        =     (U    .a2#U3#04#T1#..  .74,W1#...W4  ,W5,...W12f        (2.18) 

Conventional  finite  element  codes  utilize  three 
translational  and  two  rotational  degrees  of  freedom  at  each 
node  (each  element  has  5  4=20  d.o.f.  as  in  EPSA  ).  However 
the  masses  associated  with  rotational  degrees  of  freedom  are 
very  small,  leading  to  numerical  instability.  The  use  of 
the  aforementioned  formulation  alleviates  this  problem  since 
only  translations  are  considered  and,  in  addition,  the 
number  of  unknowns  is  reduced,  leading  to  simpler  and  more 
efficient   computations. 

It  must  be  pointed  out  that  in  order  to  model  the 
bending  behavior  at  the  edges  of  the  shell,  a  set  of  artifi- 
cial nodes  has  to  be  created.  The  finite  element  grid  will 
then  consist  of  the  nodes  defined  in  the  input  file  plus  the 
artificial  nodes  along  the  boundary  of  the  sheet 
(figure   2.3). 

2.      Strain    Disp  lacement    Relation 

The  principle  of  virtual  work  deals  only  with 
strains.  Since  the  finite  element  approximation  gives  the 
displacement  at  each  point,  the  equations  relating  the 
strains  to  the  displacements  are  needed.  In  this  study,  the 
Donnell-Vlasov  nonlinear  kinematic  equations  of  shell  theory 
are  employed,  and  the  strain-displacement  relations  are 
described   in  table   I. 

Using  equation  (2.15)  in  the  equations  detailed  in 
the  previous  section  will  give  the  finite  element  approxima- 
tion   of    strains    in    terms   of    nodal   displacements. 
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TABLE    I 
Donnell-?lasov   Shell   Equations 


9ul  _i_   ^i  JL      I     2 

ell  "  h13C1  +  h1h2  dc,2  u2  +  RL  +  2  *l 

C22  "  h29C2       hxh2   3^  ul       R2       2  y2 

3u2  3u.  ^      3h  ^^       3h2 

2e12  "  h^  +  h^  -h^^Ur'h^9^U2  +  V2 

34>,  ,       3h  3$  3h, 

k       -  r-  +  — 4>        k  — r~  +  — 0 

Kll       1^3^       hxh2  aq     2     *22       h23C2       hxh2   3£2  vl 

Zk12  "  h13^1  +  h2342  "  hxh2   di2  91  "  hxh2   3^  92 


„.  3w  3j* 

Where         <b     -  -  <t> 


1  hL3^  v2  h23^2 

„'£     are  local   coordinates 
l     2 

I      h1  ,h2    are  scaling    coefficients. 


3«      Shell  Constitutive    Relations 

The  shell  constitutive  relations  relate  the  stress 
resultant  rate  vector  to  the  shell  strain  rate  vector.  In 
matrix  terms: 
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(S)     =    [D]{3} 


(2.  19) 


Where  [ D  ]   is  the   Elasto-Plastic  tangent   stiffness    matrix. 


Figure   2.3        Modal  Points   Organization. 

The  stress-strain  relation  used  in  EPSA  differs  from 
the  classical  Elasto-Plastic  theory  in  that  the  formulation 
involves  shell  stress  resultants  rather  than  stresses  at 
points  throughout  the  thickness  of  the  shell.  In  other 
words,  EPSA  uses  relation  2.19  integrated  over  the  thickness 
of  the  shell.  Thus,  there  is  no  need  to  compute  the  stresses 
throughout  the  shell,  which  results  in  a  significant  savings 
in  storage  space  and  procsssing  time.  However,  the  stress 
resultants  Ni>and  m\  .  cf  the  shell  theory  are  not  sufficient 
to  describe  the  state  of  stress,  and  certain  higher  order 
moments  must  be  combined  with  the  strass  resultants  to  form 
the  dynamic  variables  of  the  problem.  The  coefficients  for 
the  integrated  constitutive  equation  have  been  determined 
using   experimsntal   results    [Ref.    5]. 
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These    constitutive    equations  consist    of    a    yield   condition,    a 

strain   hardening   law    and   a    flow   rule: 

-The    yield   condition   indicates    whether   part   of   the   shell   has 

started   to    yield.       (figure    2.4) 

-The    strain   hardening  law  gives   the    evolution   of    stresses   in 

the    shell  after    plasticity    is   reached. 

-The    flow      rule    gives  the   plastic      strain   rate   in      the   shell 

after  plasticity. 


W2 


"Z. 
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a  Oyield 


/* 
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Figure    2.4        Yield    Situation   in   the   Shell. 


4.      Solution   Procedure 

EPSA  uses  an  explicit  central  difference  scheme  to 
solve  the  virtual  work  and  fluid  loading  equations  detailed 
in  section  B.  As  indicated  in  appendix  B,  the  explicit  time 
integration  procedure  requires  a  small  time  step  and  is  only 
conditionally  stable.  However  when  stable,  it  always 
converges  to  the  exact  solution,  as  opposed  to  implicit 
schemes  that  are  unconditionally  stable  but  may  lead  to 
unrealistic  results.  Furthermore, ia  problems  involving  the 
treatment  of  shocks,  accuracy  requirements  preclude  the  use 
of  large  time  steps.  The  central  difference  scheme  seems, 
therefore,    particularly   optimal. 
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Assuming  we  know  the  solution  at  tine  t  ,the  central  differ- 
ence scheme  applied  to  equation  (2.17)  gives  the  solution  at 
time    t+At    : 


fc*At  t  4 

{V}.       =       {7}.     +     _At_   ({PJi    -    X.  ^Fhc  >  (2.20) 


Mi 


k=l 


Where  & .  is  the  mass  of  node  i  ,  {PJi  are  the  externally 
appied  forces  and  (F)  are  summed  over  all  the  elements  k 
framing   node  i. 

The  formulation  of  the  equations  is  in  the  initial 
configuration  and  all  equations  are  solved  in  the  initial 
geometry  in  accordance  with  the  total  Lagrangian  formulation 
(Ref.    6]. 

E.       EPSA   CAPABILITIES 

The  structure  to  be  analyzed  is  conceptually  divided 
into  constitutive  parts  called  "shaets."  Each  sheet  is  a 
curved  sectioa  of  a  shell  with  an  arbitrary  number  of  nodes 
and  elements  (figure  2.5)  Its  shapa  is  limited  to  a  surface 
that  can  be  described  by  a  smooth  continuous  function 
without  discontinuities  in  its  slope.  The  elements  within 
the  sheet  are  limited  to  a  rectangular  organisation  (figure 
2.5)  . 

Thus,  a  cylinder  with  end  caps  would  consist  of  three 
sheets:  a  cylindrical      sheet      and    a      sheet      for    each      end 

(figure    2.5).  Three   sheets     are   required      because   of      the 

slope  discontinuity  at  the  edge  between  the  cylinder  and  the 
end   caps. 

Dividing  the  structure  into  sheets  isolates  the  diffi- 
culties associated  with  the  boundaries  into  a  set  of  artifi- 
cial nodes  along  the  perimeter  of  the  sheet.  When  several 
sheets  are  required  ,EPSA  takes  care  of  the  compatibilities 
of   displacements,    rotations   and   moments   along   the    edges. 
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Figure   2-5        Sheet   Organization. 

Thus,  any  arbitrarily- shaped  structure  can  be  analyzed 
using   EPSA    by   dividing  it   into  a    number   of    sheets. 

Two  options  for  choosing  elements  are  available  in  EPSA. 
The  first  option  exists  to  employ  \  generalized  quadrilat- 
eral element.  The  second  oprion  exists  to  employ  a  rectan- 
gular element  and  uses  a  specialized  formulation  for  this 
type    of   element. 

The  coordinates  which  can  be  either  cartesian  or  cylin- 
drical always  lie  within  the  plans  of  the  sheet.  The  z 
direction  lies  in  a  positively  outward  direction  normal  to 
the  sheet.  Each  sheet  contains  its  own  local  coordinate 
system,  there  is  no  global  coordinate  system  when  multi 
sheets   are    merged    (figure  2.6). 

Prior  to  generating  a  finite  slement  mesh  for  a  sheet 
one  must  establish  the  side  numbers  of  the  sheet.  The  side 
numbering  scheme  is  arbitrary  as  to  the  choice  of  sheet 
number  one.  However  the  specification  of  sides  1  to  4  must 
proceed   in      a   counterclockwise   direction      when  the      sheet   is 
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Figure    2.6        Coordinate   System. 

viewed  from  the  positive  z  direction.  At  the  intersection 
of  side  1  and  side  4,  element  1  is  first  defined.  Then  the 
node/element  number  is  incremented  by  one  until  it  reaches 
side  2.  Then  it  returns  to  side  4  and  continues  the  count 
for   the   next  row   of    nodes/elements. 

Thanks  to  the  exclusive  use  of  quadrilateral  elements 
and  to  the  specific  numbering  procedure,  the  table  of 
connectivity  is  implicitly  defined  *hen  generating  the  nodal 
points,  thus  simplifying  considerably  the  input  require- 
ments. 

The  inclusion  of  structures  internal  to  the  cylindrical 
sheet  is  enacted  in  EPSA  through  the  use  of  internal  sheets. 
Structures  internal  to  a  cylinder  are  therefore  modelled  as 
individual  sheets  having  the  saae  capabilities  as  any 
general   EPSA  sheet. 

Two   types   of   internal   structures    are   available: 
-Sheets    (beams,      plates   or    shells)         whose   connection   to  the 
cylindrical  shell  lies  parallel   to   the   axis   of  the   cylinder. 
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-Shaets  (beams,  plates,  shells)  whose  connection  lies  in  the 
circumferential    direction   of  the   cylinder. 

Furthermore, in  order  to  modal  internal  equipment 
(machinery,  etc  )  which  does  not  deform  but  contributes  to 
the  inertia  of  the  system,  concentrated  masses  can  be  input 
at    nodal   points. 

The  user  aust  be  aware  that  the  previously  discussed  DAA 
is  only  implemented  for  cylindrical  stuctures.  Prior  to  the 
finite  element  analysis  the  user  must  compute  the  added  mass 
(virtual  mass)  matrix  defined  in  aquation  (2.4).  This  is 
done  by  using  the  ACCESION  program,  which  creates  a  virtual 
mass  file  that  EPSA  reads  when  computing  the  flu  id- structure 
interaction. 

Finally,  EPSA  in  its  latest  version  takes  local  cavita- 
tion into  account.  When  the  total  pressure  computed  by  EPSA 
is  found  to  be  negative,  it  is  simply  set  to  zero  since 
water  cannot   withstand  any    tension. 

F.       USING    EPSA 

The  input  file  for  EPSA  can  be  generated  either 
directly,  or  via  the  interactive  program  PEEPSA  that  prompts 
the  user  to  give  the  input  data  .  When  the  input  file  is 
created,    all  the  data  are   in  free   format. 

The  nodal  points  can  be  generated  semi-automatically 
(see  the  user's  manual),  and  this  option  is  very  helpful  and 
time  saving  when  generating  big  models.  EPSA  can  be  run 
interactively  for  small  models  or  on  batch  for  bigger  jobs. 
For  instance,  a  cylindrical  shell  with  1440  elements  and 
1517  nodal  points  takes  0.0  129  sec.  of  VAX  CPU  per  time  step 
per  element.  The  whole  model  would  take  about  half  an  hour 
for    100    steps. 

EPSA  creates  an  output  file  in  which  all  input  data  is 
echoed,    and   outputs    the   pressures,    stresses,    strains,    veloc- 
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ities,      displacements  at      the   nodal    points   requested      by   tha 
user    in   the   input  deck. 

The  value  of  the  time  increment  At  should  be  selected  so 
that    : 

At      <      1/2  6m.n(  P/E)"5  (2.21) 

Where<5    .is  the      smallest   disxance   between   the      midpoints   of 
mm 

opposite   sides    of  an    element,      for   all   elements   of   the    sheet 

(figure    2.5),      and    1/2  is  a    safety   factor.         In   other   words, 

tha    time   step      increment   has   to    be  less    than      the   time   taken 

by   a    wave  to  propagate   from   an   element,   to   another.      In    equa- 
ls     . 
tion    2.21,     (E/p )       is    simply    tha   wave    speed   in   the    material. 

The  virtual  mass  array  (VMA)  is  created  on  unit  20, 
therefore  one  should  not  use  this  unit  for  any  other  purpose 
than    READ  operations. 

Because  of  the  simple  organization  of  its  input  file, 
EPSA  has  been  found  fairly  easy  to  use.  The  user  can  perform 
major  changes  in  the  model  very  quickly  and  efficiently.  The 
ACESION  module  has  been  found  to  wock  well  except  for  cylin- 
ders of  large  dimensions  (L=1400",  D=240")  for  which  the 
size  of  the  virtual  mass  array  grows  unexpectedely  from  a 
reasonable    4  blocks    to    190    blocks   of    VAX   memory   size. 

However,      EPSA    has  been    designed    for   some   specific   types 
of   f luid-stucture      interaction   analysis   and      its    limitations 
should   be   pointed  out: 
-  Only   beam  type   stiffeners    can   be   considered 

The  fluid  structure  interaction  is  only  enacted  for  a 
circular  cylindrical  sheet,  which  takes  away  much  of  the 
flexibility  the    program. 
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III.     EPSA/   PATRAN3G    INTERFACE 

A.       INTRODOCTIOH    TO    COLOR   GRAPHICS    SYSTEHS 

In  dealing  with  the  response  of  a  structure  to  a 
loading,  quantities  such  as  stresses,  strains,  velocities 
and  displacements  must  be  analyzed  at  a  number  of  nodal 
points,  which  makes  the  interpretation  of  computer  outputs 
very  tedious.  Color  graphics  systems  offer  an  effective 
solution  to  this  problem  by  providing  a  global  representa- 
tion of  a  quantity  distribution  over  a  stucture  rather  than 
a  discrete  representation  given  by  a  computer  output.  A 
color  graphics  system  consists  of  an  interface  between  the 
computer,  the  user  and  the  color  terminal.  A  typical  system 
would  be  a  software  package  that  allows  the  user  to  create 
models  on  the  screen  as  well  as  to  display  any  data  in  terms 
of  color  intensity.  It  is  known  that  a  picture  is  worth 
several  hundreds  of  words.  Therefore,  merging  a  finite 
element  program  with  a  color  graphics  system  would  be  a 
major   improvement   in    engineering   analysis. 

PATRAN-G  [  Ref.  7]  is  a  color  graphics  system  specifi- 
cally designed  for  finite  elements,  it  permits  the  engineer 
and  the  computer  to  work  together  cowards  the  creation  of  a 
model.  The  designer  creates  an  image  on  the  screen  and 
PATRAN  automatically  translates  the  physical  model  into  a 
standard      finite      element   input      deck.  Another      important 

feature  of  PATRAN-G  is  its  ability  to  assist  the  user  in  the 
interpretation  of  results  through  its  post-processing  facil- 
ities which  include:  the  deformed  geometry  with  magnified 
deformations, the  color  coding  of  elements  based  upon  any 
response  parameters  such  as  displacements,  stresses  and 
strains.      In  particular,    the  contour    levels   of  the 
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aforementioned  quantities  can  be  superimposed  on  a 
3-dimensional  image  of  the  model,  thus  allowing  for  a  global 
analysis  of   a   complex  structure. 

B.       BERGIIG    EPS  A    AND    PATRAN-G 

As  described  in  chapter  I  ,  the  structures  under  study 
have  fairly  short  EPSA  input  files.  Furthermore,  in  this 
study  dealing  with  fluid-structure  interactions  on  a 
circular  cylinder,  only  structures  that  consist  of  one  sheet 
are    considered.  Therefore,      because  of  the      simplicity   of 

both    the   input    file    and  tha    stuctura    under   study,      there   was 
no   need   to    design  a    module      converting   a   PATRAN-G    model  that 
is   created   on   the  screen  into  an   EPSA    input   fila. 
The    remaining  tasks    were  the   following: 

-Display  the  original  finite  elemaat  model  defined  in  the 
EPSA    input    file    on    the  screen    (original   geometry)  . 

-Display  the  nodal  points  and  alamant  results  that  are 
computed   from  EPSA    on  the   screen    (postprocessing). 

1 •      Original  Geometry 

The  input  of  a  finite  elanent  model  into  PATRAN-G 
can  be  done  by  creating  a  "neutral"  file,  as  described  in 
PATRAN-G  teruinology.  The  neutral  file1  is  intended  to 
provide  a  simple  link  between  PATRAN-G  and  the  outside 
world.  It  is  written  entirely  in  80  character  card  images 
and  all  the  data  is  organized  in  saall  "packets"  of  two  or 
mora      card    images.  Each    packet      oontains   the      data    for     a 

fundamental  unit  of  the  model  such  as  node  coordinates  or 
elaments  definition.  Since  our  only  purpose  was  to  display 
the   original  geometry  of   the   structure,      a    limited   number   of 


1  Additional    information    about   nauzral   files    can   be    found 
in   the  PATRAN-G    user's   manual   [Ref.    7] 
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data    packets    (4)    was   to    be    created: 

-File  title  (packet  25)  :  two  cards  containing  the  user 
title. 

-Summary  data  (packet  26)  :  two  cards  containing  the  number 
of  nodes,  elements  and  the  date  and  time  at  which  the 
neutral   file  was  created. 

-Node  data  (packet  1):  contains  all  information  concerning 
nodes:  node  number  and  coordinates  in  a  global  coordinate 
frame. 

-Element  data  (packet  2)  :  contains  the  connectivity  table 
for   the   finite    element   model. 

-End    of   file      (packet  99)       :   end-of-file   cards. 

We  have  seen  in  the  presentation  of  EPSA  in  chapter  I  that 
the  nodes  are  defined  in  a  local,  sheet-attached  coordinate 
system,  that  artificial  nodes  are  created  along  the  edges  of 
the  sheet  to  aodel  the  bending  behavior,  and  that  nc  connec- 
tivity table  was  input.  Therefore,  the  translator  module 
that    would    be   created  had  to: 

-skip  the  set  of  artificial  nodes  and  properly  renumber  the 
grid 

-perform   a   change  of   coordinates   for    all   local   data 
-generate  the   connectivity    table. 

It  was  decided  to  employ  a  modular  design  in  which 
each  routine  would  perform  a  specific  task.  A  modular  design 
would  allow  further  changes  to  be  made  quickly  and  effi- 
ciently by  modifying  only  the  relevant  routines.  The  imple- 
mentation in  EPSA  was  made  using  a  series  of  "calls",  thus 
minimizing  the  risk  of  interference  with  the  finite  element 
computations. 

The  translator  module  crsated  was  made  of  four 
routines: 
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-PRELIM  :    computes  the  number    of    nodal    points,    the    number 

of   elements   and    displays   the   first   two   data   packets     (25,26) 

-SHEETF  :      scans      through   the     aodas,      skips      artificial 

nodes,   renumbars   the    grid,       performs    the   required   changes   of 
coordinates   and    displays  the   node   data    packet    (3  1) 

-SHCONN  :      scans    through    the   nodas,      connecting   each  node 

to   the   elements    it    belongs    to     and   displays   the   element   data 
packet    (02)    on    the    neutral    file. 

-ENDNEU  :    displays  the   and-of-fila    packet    (99) 


2  -      Using  the  Translator   Modula 

The  translator  calls  were  iaplemented  in  the  routine 
REPORT  of  EPSA.  Any  run  of  EPSA  creates  a  neutral  file  on 
unit  19.  The  aeutcal  file  name  is  therefore  FOR019.DAT  if  no 
"ASSIGN"  statement  has  been  issued  prior  to  the  computer 
run.  The  finite  element  model  might  than  be  displayed  on  the 
graphic  terminal  (Ramtek,  Tektronix)  via  the  neutral  input 
mode  of  EPSA  (see  [Ref.  7]  for  more  details).  \n  example  of 
finite  element  model  output  on  Tektronix  4014  is  given  on 
figure   3.1    . 

3-      Implementation  of  Postprocessing   Capabilities 

Postprocessing  capabilities  exist  to  assist  the  user 
in  the  interpretation  of  computer  results.  It  is  simply  a 
process  of  generating  displays  and  raports  based  upon  a 
combination   of    geometry  and   the  results   of   an   analysis. 

The  results  cf  analyses  ara  transmitted  to  PATRAN-G 
for  postprocessing  via  "neutral  rasults  filas"  as  decribed 
in      PATRAN-G  terminology.  Unlike   the      model  neutral      file 

described  in  the  previous  section,  results  files  are  in 
binary  rather  than    in  card    image    form. 
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Figure  3.1        Example  of    Finite   Element   Hodel   Display. 

One  can    distinguish    between   two   major   kinds   of   post- 
processing   displays:    deformed   geometries  and   element 
guantities. 

a.      Deformed   Geometry 

A  displacement  results  data  file  reguired  by 
PATRAN-G  had  to  be  created.  Again,  the  module  created  had  to 
skip  the  artificial  nodes,  perform  changes  of  coordinates  as 
well  as  to  write  the  nodal  diplioements  in  the  neutral 
results      file.  The     displacement    results      data      file      was 

created  in  module  NEUDISP.  Its  organization  is  given  on 
table   II. 

A  small  module  PLOTDISP  that  decides  at  which 
time  steps  the  results  should  be  output  was  created.  The 
"call"  to  PLOTDISP  was  implemented  in  module  COMPUTE  of 
EPSA.  The  overall  structure  of  the  translator  module  is 
presented  on  table    IV  at   the  end   of   the    chapter. 
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TABLE    II 
Organisation  of   Displacement   Results   Pile 

Column  Content 


1  X   displacement   in   global    coordinates 

2  Y   displacement    in   global   coordinates 

3  Z  "  "  " 

U  Displacement    normal   to   the   shell    without 

rigid   body   mode. 
5  Velocity   normal   to   the   shell 


A  sample      of  deformed  geometry      output  on   Tektronix      U014   is 
given   in   figure    3.2. 

b.      Element    Quantities  Postprocessing 

Any  element  related  quantity  can  be  the  target 
of  postprocessing.  In  general  thess  types  of  quantities  are 
the  results  of  finite  element  computations  supplied  to 
PATRAN-G  through  the  neutral  element  results  file.  The 
neutral  element  results  file  is  different  from  the  neutral 
displacements  results  file  detailed  in  the  previous  section, 
however  it  shares  a  similar  format  "m  Ref .  7].  Element  quan- 
tities such  as  stresses  in  local  and  global  coordinates  and 
von  Mises  stresses  are  computed  in  module  NEUSTRE  whose 
overall   structure   is    similar  to   NEUDISP    described    earlier. 
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Figure  3.2       Deformed   Geometry  Output. 

Ths      organisation   of  .   the  neutral      element      results   file     is 
given   on   table    III. 

As  described  in  chapter  I,  EPS  A  does  not  compute 
tha  stresses  through  the  thickness  of  the  shell.  Instead  one 
uses    the  integrated    quantities   of   the    shell   theory   [Bef.    8] 


^  h/2 

■±j-  jaij  dt 
Jh/2 


and 


01  j  tdt 


One  can  expect  the  strasses  on  the  shell  to  be 
maximum  at  the  extreme  fibers.  NEUSTRE  computes  the  von 
Mises  stresses  at  the  top  and  bottom  fibers  and  writes  the 
maximum  value  in  the  neutral  file.  At  the  surface  of  the 
shall  no  shear  stresses  are  involved.  Assuming  a  linear 
distribution  of  bending  stresses  and  a  unifori  distibution 
of   membrane   stresses,    one   can    easily    derive    : 
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aij     =      Nij  /h   *       $0  6/h2  <3-1) 

The  first  term  of  the  previous  equation  is  the  membrane 
force  contribution,  the  second  is  the  bending  moment  contri- 
bution. The  von  Mises  stresses  are  then  computed  using  the 
well-known    relation: 


a 
vm 


■  V*^      ~ai   °2     +°2        )  <3'2) 


In  a  similar  way  to  the  displacements  results 
file,  a  module  PLOTSTRE  that  de-ides  whether  or  not  to 
output  the  element  results  was  created  and  called  from 
COMPUTE.  The  overall  structure  of  tae  translator  is  given  on 
table   IV. 


TABLE    III 

Organisation   of  the  Heutral   Element   Results    File 

Column        Content        Description 

Element   local    stress,    direction   i 
ti  (i  ti  ti  -j 

Element  global   stress,    direction  x 

II  It  II  II  y 

Element  global   shear,    direction  xy 
von    Mises   stress 


22 

st re, 1 

23 

stre, 2 

25 

stre, 4 

26 

stre, 5 

27 

stre, 6 

31 

von 
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c.      Displaying  the    Results   Quantities   on    the   Screen. 

The  EPS  A  input  deck  has  been  modified  so  as  to 
create  the  results  files  (displacemsnts,  elements)  required 
by  PATRAN-G.  At  the  end  of  the  second  input  card  the  user 
specifies  the  number  of  displacement  results  and  the  number 
of  element  results  files  to  be  created  (at  least  one). 
Obviously,  those  two  inputs  are  also  in  free  format.  The 
neutral  results  files  (elements,  displacements)  will  be 
generated  at  equal  time  intervals  as  the  computation 
proceeds.  The  results  corresponding  to  the  last  time  step 
are    always   output. 

The  element  results  file  is  created  on  unit  16 
and  the  displacement  results  file  on  unit  18,  thus  corre- 
sponding to  files  FOR016.DAT  and  FDR018.DAT  respectively.  A 
new  version  of  those  files  is  created  each  time  an  output  is 
requested. 

Example: 

If  five  (5)  neutral  element  results  files  are  requested  on 
the  input  card  of  a  run  of  20  steps,  five  files  FOR016.DAT;1 
to  FOR016.DAT;  5  will  be  created,  corresponding  to  time  steps 
U  to    20    respectively. 

For  the  displaying  of  elament  and  nodal  points 
results,  the  user  will  refer  to  [Ref.  7]  section  20.  The 
title  of  the  run  (first  card  of  EPSA  input  deck)  will  be 
displayed  on  the  screen  along  with  the  time  at  which  the 
results   have  been   reguested. 
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TABLE    IV 
Structure  of  the   Translator 


V 


CALL    PL3TDISP 


no   output 
requested 


CALL    PLOTSIRE 

no    output 
requested 


computation 
proceeds 


•output   requested- 


■>— output   requested 

1 , 

CALL    NEOSTRE 


V 


ICALL    PRELIM 

1 

EjCALL    SHEETF 

PI 

0  ICALL    SHCONN 

Rl 

T  ICALL    ENDNEU 


sea  ns  through | 
the  elements-l 
change  coord  j 
inates-  ~| 

compute  Von-  j 
Mises  stress  | 
-display  J 

element  | 

results  | 


T 


RETURN 


-CALL    NEDDISP 

I 

scans    through 
the   grid- 
changes   coor- 
dinates-find 
node    with    max 
deformation- 
display    first 
card- 
scans    through 
the    grid- 
display   nodal 
results 


t 


RETURN 
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IV.    DESCRIPTION    OF    THE    ONDER HATER    EXPLOSION 

A.       PRESENTATION 

An  explosion  is  a  chamical  or  nuclear  reaction  in  a 
substance  (water)  that  converts  an  original  material  into 
other  products  plus  a  significant  amount  of  energy.  The 
process  of  the  explosion  produces  very  high  temperatures  and 
pressures  and  occurs  with  extreme  rapidity.  As  the  result  of 
the  explosion,  the  initial  mass  of  explosive  becomes  a  very 
hot  mass  of  gas  at  tremendous  pressures;  it  is  then  obvious 
that    these   conditions  will   affect   the    surrounding   medium. 

The  fact  that  the  water  is  CDmpressible  leads  to  the 
conclusion  that  the  pressure  applied  at  some  location  in  the 
liquid  will  propagate  through  it  as  a  wave  disturbance 
[Raf.    9].  It      must     be   pointed      out      that      the      pressures 

involved  in  an  underwater  explosion  are  usually  so  large 
that  the  wave  velocity  cannot  be  assumed  independent  of 
pressure.  Thus,  the  small  amplitude  wave  theory  detailed  in 
[Ref.  9]  does  not  apply  and  this  will  be  the  source  of  many 
complications  in   describing    the   behavior   of   the   shock   wave. 

The  first  cause  cf  disturbance  to  the  water  in  an  under- 
water explosion  is  the  occurrence  of  the  pressure  step  at 
the  water  boundary.  Immediately  upon  its  arrival,  the  pres- 
sure begins  to  be  relieved  by  an  intense  pressure  wave  and 
outward  motion  of  the  water.  For  explosives  like  TNT  or  for 
nuclear  explosions,  the  pressure  rise  can  be  considered  as  a 
discontinuous  step,  and  is  then  followed  by  a  roughly  expo- 
nential decay.  The  duration  of  the  phenomenon  is  of  the 
order  of  a  few  milliseconds  (more  for  nuclear  explosions)  . 
Once  initiated,  the  pressure  disturbance  is  propagated  radi- 
ally outward  as  a  compression  wave,  also  called  a  shock  wave 
because   of   the    steep    pressure  step   at    its    front. 
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Close  to  the  explosion,  the  velocity  of  the  wave  is  several 
times  the  "acoustic"  speed  of  the  small  amplitude  theory 
(c=5000  ft/sec)  ;  but  approaches  this  limit  very  rapidly  as 
the    wave  advances  outward. 

The  pressure  level  in  the  wave  falls  more  rapidly  with  the 
radial  distance  than  what  is  predicted  with  the  small  ampli- 
tude theory,  but  approaches  this  behavior  at  large 
distances. 

B.       BOBBLE    EFFECT 

As  a  result  of  the  explosion,  the  initial  mass  of  explo- 
sives becomes  a  hot  mass  of  gas  at  tremendous  pressures. 
After  the  principal  part  of  the  shock  wave  has  been  emitted, 
the  gas  pressure  is  considerably  decreased  but  is  still  much 
higher  than  the  equilibrium  pressure.  The  water  in  the  imme- 
diate region  of  the  sphere  or  "bubble"  of  gas  has  a  large 
outward  velocity  and  the  diameter  of  the  bubble  increases 
rapidly.  The  expansion  continues  and  the  internal  gas  pres- 
sure decreases  gradually,  but  the  motion  persists  because  of 
the  inertia  of  the  outward  flowing  water.  When  the  gas 
pressure  falls  below  the  equilibrium  value,  the  pressure 
defect  brings  the  outward  flow  to  a  stop  and  the  boundary  of 
the  bubble  begins  to  contract  at  an  increasing  rate.  The 
inward  motion  continues  until  the  oompressibility  of  the  gas 
reverses      the   motion.  Thus,        the      inertia   of      the      water 

together  with  the  elastic  properties  of  the  gas  provide  the 
necessary  conditions  for  an  oscillatory  system,  rhe  oscilla- 
tions of  the  gas  sphere  may  persist  a  number  of  cycles,  ten 
or  more  oscillations  having  been  detected  in  favorable 
cases. 
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C.       SURFACE    EFFECTS 

In  the  case  of  explosions  occurring  at  shallow  depths, 
surface  effects  will  complicate  the  aforementioned  sequence 
of  events.  When  the  shock  wave  hits  the  surface, the  atmos- 
phere cannot  supply  appreciable  resistance  by  compression. 
As  a  result,  a  reflected  wave  with  a  negative  pressure 
satisfying  the  zero-pressure  condition  at  the  surface  is 
formed  (figure  4.1).  Thus,  the  rasultant  pressure  observed 
is  the  sum  of  the  direct  and  reflected  pressures.  Therefore, 
the  reflected  wave  arriving  at  point  A  will  create  a  sudden 
drop  of  the  pressure  to  a  smaller  value.  This  is  known  as 
the  "cut-off"  phenomenon,  typical  of  free  surface  effects 
(figure   4.2). 


explosion 


bottom 


Figure   4.1        Surface  Effect  on  a  Shock  Have. 


D.       PRESSURE   DETERMINATION 

As   detailed    in    a    previous   section,    the   pressure    decay   at 
any    point   is  roughly   exponential   so   that   it   can   be    written: 


Pt  (A)      =     PQ  (A)    e 


ve 


(4.1) 
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It  has  been  found  that  the  fundamental  descriptors  cf  an 
underwater  explosion  attack  are  the  charge  size  (equivalent 
weight  of  TNT)  and  the  charge  staadoff  (shortest  distance 
between  charge  and  target)  .  Theoretical  developments  about 
the  spherical  wave  detailed  in  [ Ref .  10]  have  shown  that  the 
peak  pressure  Pm  at  any  point  can  be  reasonably  approximated 
by: 


Pm     =       K 


V3 


VJ     A 
l(W    /R? 


(*.2) 


where   W      is   the    charge     size   in   pouads  of      TNT   and   R      is   the 
standoff   distance   in    feet. 
It   has   been   shown  as    well: 


V3    V3     a2 
9       =      K2H    (W  /R) 


(4.3) 


K.rK  2wK  » A2  are  empirically  determined  factors  that  depend 
on  the  type  of  explosives  used.  Their  values  for  several 
types   of  explosives    are   given   on   table   V. 
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Figure   4.2        Cut-Off   Phenomenon 
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TABLE    7 

Explosion   Parameters 


HBX-l 

TNT 

PENT 

NUKE 

Pnax                        K 

22347.6 

22505 

24539 

4.33X106 

*1 

1.144 

1.18 

1.194 

1.13 

Decay  Constant  fU 

L. 

.056 

.058 

.052 

2.274 

h 

.247 

.185 

.257 

.  22 

E.       THE    EXPLOSIOH   IH    EPSA 

EPSA   features   two  different   way    of   describing    underwater 
explosions: 

-A  discrete  form  in  which  the  user  inputs  the  pressure 
history   at    a  finite    number    of    times. 

-A  functional  form  that  uses  equations  4.2  and  4.3  .The 
program  requests  then  the  various  coefficients  and  parame- 
ters   describing   the    explosion. 
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Figure  4.3       Incident   Pressure  Decay. 
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The  presence  of  free-surface  effects  can  also  be 
accounted  for  with  the  input  of  a  cut-off  time  after  which 
the    incident   pressure  is   set   to   zero    (figure    4.3)  . 
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V.     ANALYSIS    AND    fLISOLTS 

A-       HODEIS    STUDIED 

In  order  to  compare  the  results  of  the  numerical  anal- 
ysis with  the  experimental  data,  all  attention  was  focused 
on  the  Explosive  Power  Meter  (EPM>  model  for  which  field 
tests  had  been  conducted.  The  EPM  model  is  a  stiffened 
cylinder  with  end  caps  whose  dimensions  are  given  in  figure 
5.1    . 
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Figure   5.1        Explosive  Power   Meter. 
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Rather  than  modelling  the  end  caps  with  additional 
sheets  ,  it  was  found  to  be  mors  efficient  to  take  into 
account  the  behavior  of  the  end  caps  using  two  rigid  end 
blocks    (figure    5.1).  The   effects   of   the      explosion   on   the 

deformation  of  the  end  blocks  is  negligible  compared  to  the 
deformation  of  points  located  outside  of  the  end  blocks. 
Therefore,  it  was  assumed  that  the  notions  of  the  end  blocks 
were    pure  rigid    body    displacements. 

In  order  to  gain  some  insight  into  the  influence  exerted 
by  the  stiffeners  and  the  end  blocks,  a  preliminary  analysis 
was  conducted  on  a  ring  stiffened  cylinder  without  end 
blocks   as   well   as  on   an   unstiffened   cylinder    . 

In  addition  to  the  study  of  the  EPM  model,  the  influence 
that  the  location  of  the  stiffeners  had  on  the  deformations 
throughout  the  shell  was  evaluated.  By  performing  a  compara- 
tive analysis  of  the  deformations,  it  was  intended  to  mini- 
mize   and  control  the    damage    caused   to   the   structure. 

The  cylinders  tested  were  subjected  to  an  explosion 
occurring  at  the  distance  R=  200"  from  the  cylinder.  The 
location  of  the  explosion  was  symmetrical  with  respect  to 
the  longitudinal  and  transverse  axis  of  the  cylinder  (figure 
5.2) .  A  spherical  type  TNT  explosion  of  intensity  W  =  50  lb 
was  selected  for  the  study.  It  was  therefore  determined  by 
the  following  parameters  (chapter  17)  : 
kx   -     1.  18  A2     =   -.  185 

Kx  =    22505         K2     =      . 058 

The  symmetry  with  respect  to  the  xy  and  yz  plane  has 
been  taken  advantage  of  by  modelling  one  guarter  of  the 
model.  After  a  certain  amount  of  sensitivity  analysis  was 
performed  on  simple  grids,  a  finite  element  grid  consisting 
of  15  17  nodes  and  1440  elements  was  selected  (figure  5.3)  . 
For  each  of  the  cylinders  studied,  the  time  step  chosen  was 
egual  to     &t  =      3.  10**     s   .    The   explosion   process    was   studied 
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Figure  5.2    Explosion  Location. 


Figure  5.3   PEM   Discretization  1517  nodes,  1440  eleaents. 
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over  a  time  interval  of  800  time  steps,  that  allowed  for  a 
shock    and  after-shock  analysis. 

The  displacements  computed  by  EPSA  consist  of  a  combina- 
tion of  rigid  body  displacements  and  pare  deformations.  The 
deformation  modes  are  of  significant  interest  since  they  may 
induce  buckling  and  even  lead  to  the  destruction  of  the 
structure.  As  mentioned  earlier,  the  very  stiff  end  blocks 
have  pure  rigid  body  displacements.  The  displacement  of 
each  node  of  the  end  block  was  subtracted  from  the  displace- 
ment of  each  node  of  the  corresponding  row,  giving  the 
component  of  the  displacement    due   to    pure   deformation. 

For  each  of  the  aforementioned  cylinders,  the  deformed 
geometry  and  the  colcr-filled  contour  plots  of  von  Mises 
stresses  as  well  as  normal  displacement  were  displayed. 
Osing  identical  color  assignments,  the  progressive  gross 
evolution  of  the  previous  quantities  were  evaluated  so  as  to 
allow  for  a  comparative  analysis  of  the  evolution  of  phys- 
ical   parameters    throughout    the   shell. 

Color-filled  contour  plots  allow  for  a  global  rspresentation 
of  a  quantity  distribution  and  have  been  found  extremely 
valuable   in  the    interpretation   of   the   results. 

For  printing  and  processing  reasons,  it  was  not  possible 
to  include  color  pictures  in  this  document.  Instead,  the 
contour  plots  of  the  physical  quantities  under  study  have 
been    included. 

B.       ANALYSIS   OF    RING    STIFFENED    CTLIHDEBS    WITH    END    BLOCKS 

1.      2PM    Model 

The  contour  plots  of  von  Mises  stresses  at  time 
steps  20,  60,  100,  140,  180,  200  are  provided  on  figure  A.1 
to  A. 6  .  As  the  shock  wave  hits  the  structure,  it  appears 
as  if  the  stresses  propagate  through  the  shell  and  reach 
their   maximum  value    fairly    quickly      in   50    time  steps.      After 
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100  time  steps,  when  the  structure  is  fully  enveloped  by  the 
shock  wave,  the  region  close  to  the  end  block  becomes 
heavily  stressed  (figure  A. 4).  in  addition,  some  concentra- 
tion of  stresses  at  the  locations  of  the  stiffeners  can  be 
seen  (figure  A. 5).  At  later  time  stsps,  the  pressure  becomes 
decreased  and  there  is  a  relaxation  of  stresses.  However, 
tha  region  close  to  the  end  block  as  well  as  some  spots 
located  around  the  stiffeners  remain  heavily  stressed,  which 
may   indicate  local    buckling     (figure   A. 6) . 

2  •      Controlled   Deformations 

In  crder  to  obtain  a  more  uniform  distribution  of 
displacements,  the  stiffeners  have  been  shifted  towards  the 
end  blocks.  The  time  history  evolution  of  the  displacements 
was  studied  over  an  interval  of  803  time  steps  for  the  EPM 
model  as  well  as  for  the  model  with  shifted  stiffeners 
called  EPM2.  The  contour  plots  of  normal  displacements  at 
time  steps  200,  400,  600,  800  for  both  models  are  provided 
on   figure  A. 7   to   A. 14 

The  EPM  model  shows  a  significant  concentration  of 
deformations  occuring,  even  at  late  time  steps  (figure 
A. 10),  indicating  a  possibility  of  buckling.  Although  unex- 
pected, the  fact  that  the  region  close  to  the  end  blocks 
undergoes  the  most  severe  deformations  has  been  confirmed  by 
experimental  data.  A  possible  explanation  to  this  phenomenon 
is  that  the  inertia  of  the  cross-saction  of  the  cylinder  is 
relatively  uniform  along  the  cylinder,  except  at  the  end 
blocks  where  it  jumps  to  a  much  higher  value.  This  disconti- 
nuity in  the  inertia  results  in  concentrations  of  stresses 
that    cause   the   aforementioned    phenomenon. 

The  or oss -section  inertia  of  the  EPM2  model 
increases  more  smoothly  because  of  the  distribution  of  stif- 
feners along  its  axis.  Thus,  the  ooncentration  of  stresses 
has    a   lower  magnitude  and  the     region   close   to  the    end    block 
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suffers  less  damage  than  in  the  EPM  case.  It  can  also  be 
seen  that  the  deformations  are  mors  uniformly  distributed 
along  the  axis  of  the  cylinder.  Above  all,  the  deformations 
at  the  late  time  steps  ara  not  as  large,  indicating  -chat  the 
chances  of  buckling  are  significantly  lower  for  the  EPM2 
model    (figure   A  .  14)  . 

Therefore,  by  performing  an  optimization  of  the 
location  of  the  stiffeners,  the  designer  can  counteract  the 
concentrations  of  stresses  and  tha  buckling  phenomena  that 
occur  in  the  region  close  to  the  end  blocks.  It  is  believed 
that  controlled  deformations  can  have  a  very  significant 
influence  on  the  survivability  of  a  structure  submitted  to  a 
shock   wave. 

C.       ANALYSIS   OF    UNSTIFFENED    AMD   RING-STIFFENED   CYLINDER 

It  was  decided  to  study  the  progressive  gross  responses 
of  an  unstiffened  cylinder  as  well  as  that  of  a  ring- 
stiffened  cylinder  without  end-blocks.  Both  cylinders  have 
the  same  external  dimensions  as  the  2PM  model.  The  ring- 
stiffened  cylinder  is  similar  to  the  EPM  model  except  for 
the  fact  that  the  end-block  has  been  replaced  by  a  standard 
stiffener.  The  evolution  of  von  Mises  stresses  at  time 
steps  40,  80,  100  is  provided  in  figures  A.  15  through  A. 17 
for  the  unstiffened  cylinder.  Tha  evolution  of  von  Mises 
stresses  at  time  steps  40,  80,  100,  150  is  provided  in 
figures  A. 18  through  A. 21  for  the  ring-stiffened  cylinder 
without   end-blocks. 

For  the  unstiffened  cylinder  it  is  observed  that  the 
stresses  propagate  quickly  throughout  the  shell  and  that 
within  a  hundred  time  steps  an  instability  phenomenon  occurs 
showing  the  existence  of  local  buckling  (figure  A. 17).  At 
later  time  steps  the  buckling  spreads  over  the  entire  stiff- 
ened  shell. 
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The  evolution  of  von  Mises  stresses  in  figures  A. 18 
through  A. 21  shows  that  the  ring-stiffened  cylinder  can 
withstand  much  higher  stress  levels  than  the  unsxiffened 
shell.  This  result  was  expected  since  the  stiffeners  provide 
the  stiffness  required  to  withstand  higher  loads.  At  time 
step  100  the  unstiffened  shell  is  already  subjected  tc  local 
instability  characterized  by  a  "hari  spot"  in  the  middle  of 
the  model  (figure  A.  17).  On  the  other  hand,  the  stresses  in 
the  stiffened  cylinder  are  much  more  evenly  distributed 
throughout  the  model,  with  high  amounts  of  stresses  concen- 
trated  around   the  locations    of  the  stiffeners    (figure   A. 21). 

It  can  also  be  observed  that  a  significant  concentration 
of  stresses  occurs  at  the  extremitias  of  the  stiffened  shell 
(figure      A.  21).  Recalling     that   the      end-block      has      been 

replaced  by  a  standard  stiffener,  the  cross-section  of  the 
shell  has  a  greater  inertia  at  the  extremities  due  tc  the 
face  that  the  two  stiffeners  locarad  at  the  extremities  are 
close  to  each  other.  Therefore  this  phenomenon  is  similar  to 
the  one  encountered  when  studying  the  EPM  model.  However 
the  concentration  of  stresses  for  the  ring-stiffened  cylin- 
drical shell  is  less  significant  than  for  the  EPM  model,  due 
to   a    smaller  discontinuity   in   cross-section   inertia. 
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VI-     CONCLOSigN 

A  FORTRAN  module  that  merges  the  finite  element  coda 
EPSA  with  the  colcz  graphics  system  PATRAN-G  has  been 
designed  and  succasfully  completed.  The  non-linear  elasto- 
plastic  responses  of  various  types  of  submerged  cylindrical 
shells  have   been    evaluated    using   the    EPSA/PATRAN-G   system. 

The  analysis  of  the  progressive  gross  responses  of  a 
ring-stiffened  cylindrical  shell  with  end-blocks  is  believed 
to  provide  useful  information  regarding  the  behavior  of  a 
submerged  structure  subjected  to  an  underwater  explosion. 
The  influence  of  the  location  of  the  stiff eners  on  the 
deformations  has  been  studied  and  is  also  believed  to  be  of 
significant  help  in  the  determination  of  an  optimal  design 
that  will  minimize  the  damage  due  to  an  underwater 
explosion. 

The  utilization  of  the  color  graphics  system  in  the 
interpretation  of  the  results  of  analysis  has  been  found  to 
be  an  extremely  valuable  tool.  It  is  the  author's  belief 
that  the  use  of  color  graphics  systems  will  become  increas- 
ingly important  in  the  analysis  of  oomplex  phenomena  such  as 
underwater    explosions  on   submerged   structures. 
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Figure  A. 1   Von  Hises  stresses,  time  step  20  ,  EPM  bo del. 
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Figure  A. 2   Von  Hises  stresses,  time  step  60  ,  EPH  model. 
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Figure  A. 3   Von  Mises  stresses,  time  step  100,  EPH  »odel. 
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Figure   A- 4        7on  Mises  stresses,    time  step    140,    EPfl  model. 
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Figure  A. 5   Von  Mises  stresses,  tine  step  180,  EPM  model. 
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Figure  A. 6   ?on  Hises  stresses,  time  step  200,  EPH  model. 
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Figure  A. 7        Hor«al  Displace nentsr    step   200 ,   EPH   nodel. 
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Figure  A. 8    Noraal  Displacements,  step  400,  EPH  model. 
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Figure  A.  9    Horial  Displaceients,  step  600,  EPS  model. 
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Figure   \. 10        Normal  Displacements,    step   800,    EPH   model. 
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Figure   A. 11        Hor«al  Displacements,    step   200,    EPM2   aodel. 
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Figure  A. 12        Nor«al  Displacements,    step  400,   EPH2   nodel. 
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Figure  A. 13   lormal  Displacements,  step  600,  EPM2  aodel, 
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Figure  A. 14   Hormal  Displacenents,  step  800,  EPH2  model, 
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Figure  A.  15       Von  Hises   stresses,      step   40    ,    unstif-    shell. 
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Figure  A- 16        Von  Mises  stresses,      step   80    ,    unstif.    shell. 
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Figure  A. 17       Von  Hises   stresses,      step    100,    anstif.    shell. 
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Figure  A. 18   Voa  Sises  stresses,  step  40  ,  stiff,  shell. 
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Figure  1.19   Toa  Hises  stresses,  step  80  ,  stiff,  shell. 
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Figure  A.  20   Von  Mises  stresses,  step  100,  stiff,  shell 
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Figure   A. 21        Von    Mises   stresses,   step    150,    stiff,    shell. 
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APPENDIX    § 
REVIEW    OF   NONLINEAR    FINITE    ELEMENTS 

A.  INTRODUCTION 

This  appendix  is  intended  to  give  the  reader  seme 
insights  into  nonlinear  finite  elements.  The  reader  is 
assumed  to  have  some  previous  knowledge  of  finite  clement 
theory.  The  basic  principles  of  the  theory  will  be  quickly 
reviewed,  but  the  study  will  focus  on  the  problems  that 
occur      when   daaling      with   nonlinear      theory.  Most   of      tha 

information  has  been  taken  from  [ Ref .  6]  as  well  as  from 
the    course   the    author  had  at  M.I.T.    with   K.J.    Bathe   in    1982. 

B.  THE    NEED  FOR    A    NEW   THEORY 

Considering  a  coordinate  frame  defined  by  (i,j,k),  a 
body  of  volume  V  in  which  point  i(x1,x2,x3)  is  subjected  to 
the  displacements  (u1,u2  , u3)  ,  corresponding  to  a  strain 
vector      {e}     (figure    B.1). 

In  the  following  sections,  the  superscripts  0  and  t  will 
refer  to  the  body  at  time  0  aad  t  respectively,  the 
subscripts  0  and  t  will  refer  to  tha  configuration  at  time  0 
and      t   respectively.  This  chaptar,         for      tha    purpose     of 

simplicity,  will  first  deal  with  the  static  nonlinear  anal- 
ysis   of   the   material. 

In   the    linear      theory  of   finite   elements,        one  uses   the 

well    known      Cauchy    stress   tensor      Tj_j        associated      with  the 

engineering  strain   tensor  e^  ^f^  ♦    3Uj"|  . 

2  (_3Xj         5cJ  J 

Then,    the   principle    of  virtual   work   is   written: 

where   fcR   reprasents      the   virtual   work    of      externally   applied 

forces  and  5e  is  a    virtual    strain. 
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Figure  B.1   Geoietric  Conventions, 
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(B.1) 


Equation  (B.1)  is  then  discretized  over  the  body  and  becomes 
a  set  of  integrations  over  each  of  the  finite  elaments.  In 
the  case  of  a  large  displacement,  the  volume  of  the  body 
over  which  the  integration  is  performed  might  have  signifi- 
cantly changed.  Also  notice  that  aquation  (B.1)  is  written 
in  the  original  coordinate  frame  (defined  at  t=3)  and  that 
tv  and  e^  refer  to  the  current  configuration  of  the  body. 
The  Cauchy  stresses  at  time  t*  a  t  cannot  be  obtained  by 
adding  an  increment  due  to  the  straining  of  the  material  to 
the  stresses  at  time  t  .  The  rigid  body  rotation  of  the 
material  has  to  be  taken  into  account  since  the  Cauchy 
stresses  vary  under  rigid  body  motions.  Therefore,  we  must 
perform  the  integration  of  equation  (B.  1)  over  the  unknown 
current  volume  with  respect  to  tha  original  geometry  that 
could   be   significantly   different    from   the   current    one. 

The  above  discussion  emphasizes  the  need  for  a  new  set 
of  stress  and  strain  tensors  that  would  alleviate  the  afore- 
mentioned problems  and  enable  the  integration  of  the  prin- 
ciple  of   virtual  work  to   be   performed. 
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C.       DEFINING   NEW    STRESS    AND    STRAIN    TENSORS 

1 .      Green-Lagrange   Strain   Tensor 

The   structure     introduced   earlier   is      considered  and 
two    points   A  and  B   at  t  =  0,        of  coordinates0*  (°xj       ,       °B(°x) 
are    defined   .      At   time   t,      the      body    has   been   deformed    and  A 
and    B  have   moved   to    tA(tx)     and      B(   x)     (figure   B.2). 


xl 


Figure  B.2        Displacements   Conventions- 

A  Taylor   expansion    is  used    to    express      the   coordinates    of      B 
as   a    function   of   the    coordinates   cf   A. 


or    with      d  txi=    txi   -     x±      and     d  °x±  -    °xi    -    \    ,    g 


(B.2) 


lves    : 


(d  V)  =   (aj^  )  (d  \) 


3   x. 


(B.3) 


{d^}      =     [  tx]   {dox} 


(3.4) 


where      [  £x  ]  =    (  ^i )       ,  (d1^}    =    (i  4i)     ,  {d  °x»     =    (d0Xi) 

3~xj~ 

In   other   words,      equation  (B.4)       expresses  how   a    small    fiber 

defined      by     the   vector  (d    °x}      at      t=0  has     rotated      and 
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extended   between   time  0    and    time    t      when   it    becomes      {dfcx}    . 

The    matrix    [qX]    is    called  the    "deformation   gradient"    . 

t 
The    new   length      ds    of  the   fiber   will    be    : 

(tdS)2      =        (dtx}T    {dtx}  (      T   refers      to   the      transposed 

matrix) 

and   therefore 
(^S)*      =       (C5x]fd°x})T    (C^X]{d0x}) 

(tdS)2       =        (d°x)     [qC]    {d°x}  (B.5) 

t  t        T        t 

Cocl  =  Coxl  foxl  is  called  the  "  Cauchy-Gr  een  deformation 
tensor". 

Notice  that  if  the  Cauchy-Green  deformation  tensor  is  iden- 
tity, equation  (B.5)  indicates  that  the  length  of  any  fiber 
will  not  vary.  In  other  words,  whenever  rigid  body  motions 
are  considered,  the  Cauchy  Green  deformation  tensor  is  iden- 
tity   since   the    fibers  do  not  stretch. 

The  principles  of  the  finite  element  method  will  now 
be  quickly  recalled.  Assume  that  the  displacement  (u^ )  of 
any  point  of  the  body  can  be  written  as  a  combination  of  the 
displacements   of  N    selected    points  sailed   "nodes"    : 

(V      =     L     hk(ui,k  (B*6) 

k=l 

Where  the  h  are  interpolation  functions  that  depend  only  on 
the  geometry  of  the  body.  In  addition,  the  nodes  are  chosen 
so  as  to  get  a  division  into  quadrilateral  elements  and  it 
is  assumed  that  the  displacement  of  any  point  is  only  a 
function  of  the  displacements  of  the  corner  nodes  of  the 
element    it   belongs    to.    Then    equation     (B.6)    becomes 

4 
(Ui)       =  £    h^.  (Ui)k  (B.7) 

k=l 

Recalling  that    (m±)       is   simply       (fcXi)       -    (°xi)      gives: 
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(  \\      =     (°x)       ♦  £    hk(ai)k  (B.8) 

k=l 

The    deformation      matrix   [      X]   can      be   expressed      very   simply 
using   the   previous   equation    : 

[JjX]      =    o\)     =     6±j  ♦    ^     ah3c<ilL  J6  (B.9) 


Define   now   the   "Green  Lagrange"   strain   tensor   by    : 

CSG1         =    V2       ([qC]   -    [I])  (B.  10) 

Hhere  [I]  is  the  identity  matrix. 

From  the  previous  derivations,  it  can  be  observed  that  the 
Green  Lagrange  strain  tensor  is  3  for  rigid  body  morions. 
The  Green-Lagrange  strain  tensor  rafers  to  the  body  at  time 
t  with  respect  to  the  initial  configuration.  This  is  why  it 
will    be    so    useful   in    dealing   with   large   deformations. 

Recall  that  the  ultimate  goal  is  to  apply  the  prin- 
ciple of  virtual  work  to  the  structure  under  study.  In 
particular,  having  defined  a  new  strain  tensor,  the  relation 
giving  the  virtual  Green-Lagrange  strain  tensor  corre- 
sponding to  a  virtual  displacement  (6  u^  must  be  known.  In 
the  case  of  a  linear  problem,  the  virtual  engineering  strain 
tensor   would  be: 

&£L.    =    V2     (       3ui*       3uj) 
3  3  x  j  51 k. 

In  the  case  of  large  displacements,  the  Green-Lagrange 
strain  tensor  should  be  used.  It  is  shown  in  [ Ref .  6]  that 
the  virtual  Sreen  Lagrange  strain  tensor  corresponding  to 
virtual   engineering    strains  is: 

5  oeij        s  Hi      3Si      6tem  (B.  11) 

3  xi       3~x~j 


or: 
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5  t«Lj     =         9li       °*j    So^ij  (B.  12) 

3xm       3xn 

Having  defined  a  strain  tensor  which  is  invariant  under 
rigid  body  motions,  a  stress  tensor  corresponding  to  the 
Green-Lagrange    strain  tensor  needs  to   be   defined. 

2 .      Stress    Measures 

Starting   with     the    Cauchy      stress   tensor  ,      the 

Piola-Kirchoff    stress  tensor  is   defined    : 

Osij   -£     t^in      Tmntxjn  (B.13) 


Where    °p  ;  fcp      are  the    densities   of   the    material  at    time    0   and 
t  respectively    and         (    txi^       -   C  tx  ]l      is      the   inverse   of   the 
deformation  tensor    defined    previously. 
Eguation   B. 13   can   be    easily    rewritten    in   equivalent   form: 

\m    s$*      oXimSSij&Cjn  (B  -  14 ) 

It   can   also   be    shown    by     using   the   principle   of   mass   conser- 
vation  that   : 


Op      =    tp    det[  Sx]  (B.  15) 


D.       PRINCIPLE   OF    VIRTUAL   WORK 

The  principle  of  virtual  work  of  equation  (B.1)  is 
rewritten  using  the  strain  and  stress  tensors  defined  previ- 
ously in   equations    (B.12)    and    (B.14)     : 


tR       .    /*■£    *X         S^kO3**       gx^X^ij       dV  (B.16) 


or: 
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/  02        oSijSoUjdV  (B.17) 


t_  /  tg        t_     _t 

Op 
V 

and    using  equation    (B.15) 

t        t 


Hi      =        /  0Sij60£ijdV  (B.  18) 


Thus,  the  principle  of  virtual  work  has  been  simply 
expressed  in  terms  of  a  new  set  of  stress  and  strain 
tensors,    integrated   over  the   original    volume   °V    . 

E.       THE    INCREHENTAL    CONTINUUM    MECHANICS    EQUATIONS 

In  this  section,  the  principle  of  virtual  work  will  be 
applied  to  the  structure  and  the  incremental  formulation 
using  the  Piola-Kirchof f  and  Green-Lagrange  tensors  will  be 
developed  .  Non  linear  terms  will  arise  from  the  rather 
complicated  definition  of  the  strain  tensor,  but  it  must  be 
pointed  out  that  the  new  formulation  provide  the  means  for 
the    modelling   of   large  deformations. 

Assume  that  the  configuration  of  the  body  at  time  t  is 
known,  the  configuration  at  time  t+At  must  be  determined. 
Writing    the   principle  of   virtual    work   at   time   t+At    gives    : 

=    Mti  (B.  19) 

t+*H      :    external   work  at   time   t+At 
5  o^ij    virtual    increment  in  G.L.    strain 
o^ij   :    stress    state   at   time  t+At 

Separating  between  the  known  terms  that  refer  to  the  config- 
uration at  time  t  and  the  unknown  terms  that  are  the  incre- 
ments of  stress  and  strain  between  time  t  and  time  t+At 
gives    : 
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t+At  t  «At  t 

0S  =     0  Sy  ♦   0Sij       and         e       =  oeij  +  o£ij     e  (3.20) 

oSij  ando-ij  are      simply   the    increments  in      stress    and   strain 

respectively. 

It        is      derived        in     [ Ref .    6]        that        the      increment        in 

Green-Lagrange    strain  oeij    is  made   of    a   linear   part      oeijand  a 

nonlinear  oneji^.  .       The   term      linear    refers   to  the   increment 

in    displacement      u^    . 

0£ij      =       G%.j  ♦oriij       cr       5    O^ij    =     5   (       oeg    +onij       )  (B.21) 

Using   equation    (B.21)    in  equation    (B.  19)    gives      : 

(  oSij  +    &LJ  (  *  (oeij  ♦olij  ))     3V      =     *A1R  (B.22) 

j 
Again,      separating    between    the    known    and    unknown    terms    gives 


1: 


JL    oSrfVijdv     «•[   os5o^ij       =    t+^  -JaSsiJ6oeij  (3.23) 


For      small      increments  in      displacements,        equation       (B.23) 
written   at   time    t   indicates    that     5  oeij  =      Soey   .    This   signi- 
fies   that    in  equation    (B.21)       the      aon   linear   term    is    negli- 
gible 

Then  the     constitutive    law   of     the   material      detailed   in 
chapter   II    allows   to   relate    stresses    and   strains: 

-0Sij    =    O^jrsOSij      ^    oCijrs      Oeij  <B*  2^ 

and    B.23   becomes: 

(«0%j  (fijs&s^      ♦(S%50r^dV=t^R    -    f^      #ijav  (B.  25) 


83 


The      right    hand      side     terms   of     ths      previous   equation      are 

known: 

^e^  is  the    linear   increment    in   virtaal   strain   involving  only 

known    terms 

Bsij    is   known  from   the  previous   time   step 
i>At     . 
R      is   the    work    of    external    prescribed   virtual   forces. 

The      left   hand      side   of     B.  25    is      unknown   since        Q^jand^^j 

involve   the   displacement   from   time   t    to   time   t+At    . 

F.       FIHITE    ELEMENT    DISCRETIZATION 

Equation  (B.25)  will  be  discretized  over  the  structure, 
using  the  finite  element  approximation  defined  in  section 
B.9    .  Let  N      be    the     number   of      aodes,      the     principle   of 

virtual     work      will    be      invoked      N      times,      setting      a      unit 
displacement  at    each    node  in  turn: 

5  u  k  =    1      ,        Su.  =    0      ,k*j 

A  system  of  equations  whose    unknowns    are  the   nodal   displace- 
ments '   is      obtained      Let      { A  u}      be      the  vector      of      unknown 
displacements,         {F}    be     the   vector      of  nodal   point      forces 
equivalent   to   the   internal    stresses. 
Then    equation    (B.25)    can  be    rewritten   in  matrix   form    : 

[X.]{Au}    +    C%L]{^}      =     |+/%}       -    f*F}  (B.26) 


[  QK  L]  and  [  qK^]  are  known  from  the  material  characteristics 
at  time  t  and  correspond  respectively  to  a  linear  and  nonli- 
near contribution.  It  is  therefore  possible  to  solve  equa- 
tion (B.25)  for  {Au}  .  However,  because  of  the  assumption  in 
equation  (B.24)  ,  the  exact  solution  might  not  be  reached 
immediately.  Furthermore,  depending  on  the  time  step  size, 
the    solution  process    might    even      be    unstable!        In    any    case, 
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an   iteration  for  solving   aquation      (B.26)       mast    be   performed 
until   rhe  exact    solution   is    reached    . 

A     widely   used      scheme  is      the      "modified  Newcoa      iteration" 
defined    by   the    following   aquation    and   boundary  conditions    : 

(    CSkL]  *   CSknl]  )  {Au)i.     =      T^}      -      f*B   f1  (B.27) 


and  f^u}12    I   u}1-V    {Auf      ,    with   the   initial  conditi 


ons 


TV    =    {'«»  and       f F}°     =    e?) 


l&\i}x  is   the  vector      of  incremental  aodal    point  displacements 

at   iteration  i. 

fc-At 

{  R}    is   the   vector    of     applied    loads    (constant   in    the   itera- 
tion) 
ti-Ati-l 
{  F}    is      the  vector    of  nodal    point    forces      equivalent   to   the 

stresses  at   time   t*At,    iteration    i-1    . 

At  the  first  iteration,  equation  (3.27)  reduces  to  equa- 
tion (B.26)  giving  an  increment  of  displacement  (Au}  .  Then 
a  better  approximation  of    oers  *  ^(Tii        ^s   obtained    .  gS^  is 

f  At. 

updated      to   the      new    state      of      strassas    and      becomes   o^Lj 
Equation    (B.27)       is    then   used   to    determine   the  new    increment 

2  ... 

in    displacement      {Au}      ,      and      so   on    until    the      increment   in 

tfrA  tn£t        i 

displacement  is   small  enough,        so   that         {   3}    =       (      F   }      in 

equation    (B.27) 


G.       INCLOSIOH    OF    DYNAMIC    FORCES 

If  the  loads  are  applied  rapidly,  inertia  forces  need  to 
be  considered  and  a  truly  dynamic  problem  has  to  be  solved. 
Using  d^lembart's  principle,  the  eLament  inertia  forces  are 
simply  included  as  part  of  the  body  forces.  Let  {u}  be  the 
vector  of  nodal  accelerations  and  [  t!  ]  be  the  mass  matrix  of 
tha  system.  Then  the  principle  of  virtual  work  is  written  in 
the    following  way   : 
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(C0iy{AU}    ♦      CqK^CAu   }       =         B      -CM}{AU}     -P  (B.28) 

Equation  (B.28)  represents  a  systam  of  differsntial  equa- 
tions of  second  order.  If  the  son  linear  term  C*KNL  ]  were 
negligible,  the  solution  could  be  obtained  by  standard 
procedures  foe  solving  differential  equations  with  constant 
coefficients.  However,  the  procedurss  proposed  for  the  solu- 
tion of  general  systems  of  differential  equations  can  become 
very  expensive  if  the  order  of  the  matrices  is  large. 
Therefore,  whenever  the  system  is  linear  or  nonlinear,  seme 
effective  methods  for  solving  the  equations  governing  the 
equilibrium   are    required. 

1  -      Direct    Integration    Methods 

The  essence  of  direct  integration  methods  is  based 
on  two  ideas.  First,  it  is  aimed  to  satisfy  B.28  only  at 
certain  time  intervals  apart.  Second,  a  variation  of  accel- 
eration velocities  and  displacements  is  assumed  within  each 
time  step.  The  form  of  the  assumption  determines  the  accu- 
racy,   stability    and    cost   of    each    method. 

In  the  following,  assume  that  the  initial  conditions 
(displacements,  accelerations,  velocities)  at  time  0, 
denoted  (  u  ,  u  ,  u  )  are  known.  la  the  solution,  the  time 
span  under  consideration,  T  ,  is  subdivided  into  n  equal 
time  intervals  At.  Assuming  that  the  solution  is  known  at 
time  t  ,  the  methods  of  getting  the  solution  at  time  t+At 
will    be    investigated. 

2 .      Central    Difference    Method 

In  the  Central  Difference  method,  a  finite  differ- 
ence   approximation    will   give   the    acoeleration   at    time   t  : 

u      =     1     ( *u    -    2  u   *      u   )  (B.29) 
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Writing     the      principle     of      virtual      work      at      time      t      and 
substituting  into  equation    (B.30)    gives    then: 

[*]{*«)     ♦   [OK    ]{tu}       =    {tR}  (B.  30) 


±&  t  t  t  frfc 

iCa]{    U}       =    {    R}     -     ([oK]    -    2[!1]/At2  )   {   U}    4[H]{    u}        (B.31) 

At"  A? 


The  previous  equation  gives  the  deformation  at  time  t+At 
from    the   characteristics   of    the   system   at   time   t    . 

When  [M]  is  diagonal,  which  is  frequently  the  case 
for  mass  matrices,  the  solution  at  time  t+  A t  does  not 
involve  any  triangular  factorization  of  the  matrix  [M]  , 
thus    leading  to    more   efficient   compa tations. 

The  shortcoming   in    the    use   of    the   central    difference 

method  lies   in   the   time   step   restriction:    for   stability,    the 

time    step  size      t   must   be      smaller   than    a   critical   time  step 

/£crwhich   is   equal  to   T  /u    ,    where   T      is      the   smallest   period 

cf   the  finite  element  system. 

The  central  difference  scheme  is  fairly  easy  to 
implement  for  the  integration  of  a  system  of  nor.  linear 
differential  equations.  However,  because  of  the  limitations 
of  the  time  step  ,  it  might  not  be  suitable  for  cases  when 
loads   are  varying  at   a  slow    pace. 

3 •      Implicit   Integration   Schemes 

Since  the  interest  of  this  study  lies  in  central 
difference  schemes,  this  section  will  be  limited  to  a  short 
description  of  the  fundamentals  of  implicit  time  integration 
schemes. 
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Implicit  time  inregration  schemes  use  the  principle 
of  virtual  work  written  at  time  t+  A  t  and  not  at  time  t  as 
for    the    central    difference    method    : 

ts-At  t4&t    tuAt  is- At 

[M]     {if}     ♦    [qK]    {u}  =       {    R}  (3.32) 

Again,      using   a      finite   difference   approximation   of      {u}    and 

tyAt 
replacing   into    equation      (3.32)       enables   to   solve      for    {   u}  . 

t-tAt 
Since   the   formulation   involves   the      rigidity   matrix   fo  K  ]   and 

tf-At 
the    external  work      {    R}    which   are   both      unknown,       the   system 

has    to      be   solved   in     a    similar      way    to   the      Modified   Newton 

iteration   that    was    detailed   previously. 

Implicit  time  integration  schemes  are  stable  regard- 
less of  the  size  of  the  time  step  used.  However,  if  the  time 
step  size  is  too  large,  significant  errors  can  be  accumu- 
lated  at   each   time   step,    leading   to   unrealistic   results. 

The  reader  will  find  mora  details  on  the  various 
implicit  method  in  [ Ref .  6].  let,  it  can  pointed  out  that 
implicit  methods  are  more  tedious  to  implement.  On  the  other 
hand,  a  larger  time  step  can  be  used  in  the  solution  proce- 
dure, which  can  be  of  extreme  importance  when  studying 
phenomena  over   a   significant   period   of   time. 
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APPENDIX    C 
HOW    TO    OSE    THE    TRANSLATOR    EPS A-PATR AN-G 

This  appendix  is  intended  to  explain  in  detail  the  use 
of  EPSA/PATRAN-G  post-processing  facilities.  It  is  divided 
in  three  sections  that  will  deal  successively  with  1)  the 
displaying  of  the  original  model  ;  2) the  deformed  geometry; 
3)    the   contour    plots   of   element   and   nodal    points    quantities. 

A.       DISPLAYING    THE    ORIGINAL    MODEL 

When  making  an  initial  EPSA  analysis  on  a  particular 
structure,  the  geometry  of  the  model  has  to  ba  input  into 
PATRAN-G.  As  explained  in  chapter  III,  all  the  geometrical 
information  is  contained  in  a  file  FOR019.DAT  that  is 
created  each  time  an  EPSA  run  is  lade.  The  input  of  the 
original  geometry  must  be  made  via  the  neutral  input  mode  of 
EPSA.  The  procedure,  starting  from  the  "logon"  to  PATRAN-G 
is   the   following   : 

-  Select   the   30    option 

-  Select  the  new   data   file    option    (option    1) 

-  Select   the  neutral   system     (option   4) 

-  Select   the  input    mode    (option   4) 

-  Input   the  neutal    file    name    :    FOR019.DAT 

The    original  geometry  will    then   be  displayed   on   the   screen 

It  is  often  found  convenient  to  have  a  perpective  view 
of   the   model  under    study.      In   that   oase,    the   user    should    : 

-  Issue  the   VIEW   command 

-  Select  the  rotation  about    the   absolute   axes    (option    1) 

-  Input   an   angle   of    rotation 

(23,-34,0      will   give  a  very  nice      view    but    any   angle   can  be 
input) 
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-   Issue      a    PLOT    command    to      have    tha    model    displayed      in   the 
ne  w   a  xe  s . 

An   example    of   the   procedure    is   provided   on    table    VI. 


TABLE    71 

■■ 

Finite    Eleaent    Model    Input   Procedure 

MODE?      1. GEOMETRY  MODEL  2. ANALYSIS  MODEL   3. DISPLAY  4. NEUTRAL 

SYS. 

5 

.END 

>4 

NEUTRAL   FILE?      1. CREATE  OUTPUT     3. INPUT  MODEL     3. POST-PROCESSING 

4 

,END 

>2 

INPUT   NEUTRAL  FILE   NAME 

>FOR019.DAT 

DO   YOU   UISH  TO   OFFSET   ANY   NEUTRAL   INPUT   IDS?    (Y/N) 

>N 

EPM   200   STEPS   ,N0   STIFFENERS   U*30. 

SHALL  UE   PROCEED   UITH  THE   READING  OF   THIS   FILE'?    (Y/N) 
> 

The  PLOT  command  can  be  issued  anytime  to  display  the  orig- 
inal   geometry  on   the   screen. 

When  studying  complex  models,  one  does  not  want  the 
elament  and  node  numbers  to  be  printed  along  with  the  geom- 
etry of  the  structure.  The  command  SET,  LABS  3,  OFF  followed 
by  a  PLOT  will  display  the  original  geometry  without  any 
labels   printed. 

When  a  model  has  been  input  into  PATRAN-G,  a  data  file 
PATRAN.DAT  is  created  on  the  user's  directory.  When 
connecting   with      PATRAN-G   at        a   later      time,      the      user   can 
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select  the  option  "last  data  file  "  (option3)  to  have  the 
original  geometry  displayed  on  the  screen  without  having  to 
input  the  model  again. 

B.  DEFORHED   GEOHETHY 

On  the  second  card  of  the  PATRAN-G  input  deck,  the  user 
specifies  the  number  of  PATRAN-G  displacement  data  file  and 
the  number  of  element  results  data  file.  Two  non-zero  inte- 
gers in  free  format  must  be  placed  at  the  end  of  the  second 
card  (section  II  in  the  user's  manual)  requesting  the  number 
of   displacements  and   of   elements    files   respectively. 

Assuming  that   the  user    has    made      a    200    steps    run   with    10 
output   requests      for    PATRAN-G    displacement    files,         ten    (10) 
files   FOR018.DAT   will     then    be   created   at      equal   time    inter- 
vals.        The   deformed    geometry   corresponding   to   time   step    100 
will    therefore    be  contained    in    FOR01 8. DAT; 5. 

To  display  the  deformed  geometry  corresponding  to  time 
step    100,   the  user    should   issue   the    following   commands: 

-  RUN,DEF  :    requests   deformed   geometry   option 

-  Input    the  name  of    the   displacements   file    :FOR0  18.DAT    5 

-  Select   the  PLOT  option    (option    3) 

Select      the      undefcrmed    geometry      (2)         followed      by     the 
deformed      geometry    (3).  An   example      of      the   procedure     is 

provided  on   table  711.  The   undeformed   geometry    superposed 

with  the  deformed  geometry  will  taen  be  displayed  on  the 
screen. 

C.  POST-PROCESS IHG    OF    ANALYIS    RESULTS. 

Element-related  quantities  like  von  Mises  stresses  are 
contained  in  FOR016.DAT  files,  nodal  point  quantities  are 
stored  in  FOR018.DAT  files.  As  described  in  chapter  III, 
each      column  of      those  files     contains      a   specific      quantity 
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TABLE    ¥11 
Deforned  Geometry   Procedure 

I  > 

I  nODE?   1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3. DISPLAY  4. NEUTRAL  SYS.  5. END 
I   >RUN,C0N,C0L,4 

I   CURRENT  FILE  FOR  NODAL  RESULTS  IS   FOR018.DAT 
NEUTRAL  RESULTS  FILE'?  l.NEU  FILE  2. CURRENT  FILE 

M 
|  INPUT  THE  RESULTS  FILE  NAME: 

>FOR018.DAT;5 

I   DATA  UIDTH  »     5 

FILE  TITLE  *EPM  800  STEPS 

7.5000129E-04 


DATA  VALUES  RANGE  FROM  -0.S2SE+00  TO  0.134E+00 
ASSIGNMENT?   l.AUTO  2. MANUAL  3. SEMI-AUTO  4. USE  CURRENT  LEUELS  5. END 
>4 

PREVIOUS  CONTOUR  LEVELS  USED. 
MODE?   1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3. DISPLAY  4. NEUTRAL  SYS.  5. END 
>RUN,HI,C 


(i.e.  column  31  of  FOR0  16.DAT  contains  the  von  Mises 
stre-sses).  The  reader  will  refer  to  chapter  III  for  the 
detailed  organization  of   those   files. 

Again,  assume  a  200  time  steps  analysis,  with  10  output 
requests  for  element  results  files.  The  user  might  want  to 
display  the  contour  plots  of  von  Mises  stresses  at  time  step 
100    .      In   this    case   the   following   commands   should   be   input   : 

-  RON,  CON,  COL,  31:  tells  PATR&N-G  to  look  at  column  31 
that    contains   the   von  Mises    stresses. 

-  Input   the   file    name        FOR016.DAT   5 

-  PATRAN-G  will  then  ask  for  a  color  assignment  (automatic, 
manual,  semi-automatic,  current  levels  used)  that  the  user 
will    select  according  to  his  needs. 

The    contour   plots   are  then    ready    to    be   displayed    : 

-  RUN,    HI,    FR   would    display    the   color-filled    contour    plots 

-  RUN,   HI,    CON    would    display  the   contour    plots    (color   lines) 
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The  von  Mises  stresses  are  the  most  useful  element:  quan- 
tities to  be  displayed,  but  other  element-related  quantities 
detailed  in  chapter  III  like  the  x  and  y  stresses  in  local 
or  global  coordinates  could  be  displayed  as  well  by  looking 
at  their  corresponding  column  in  the  element  results  data 
file. 

Dealing  with  nodal  point  quantities,  the  displacement 
normal  and  the  velocity  normal  to  the  shell  are  very  mean- 
ingful quantities  in  an  analysis.  They  are  respectively 
stored  in  column  4  and  5  of  the  displacement  results  files 
FOR018.DAT    . 

A  contour  plot  of  the  normal  displacement  at  time  step  100 
would   then    be   obtained  via    the   following   commands: 

-  RUN,    CON,    COL,  4  (look    at   column    4) 

-  FOR018.DAT  5  (name   of    the    file) 

-  Color   assignment    chcsen 

-  RUN,    HI,    C         or         RON,    HI,    FR 

Notice  When  the  fluid-structure  interaction  is  ON,  the 
normal  displacement  contained  in  column  4  corresponds  to 
£ure  deformations,  the  rigid  body  contribution  having  been 
taken  out. 

All  the  element  results  processing  is  implemented  in  the 
routine  NEUSTRE,  all  the  nodal  points  processing  is  imple- 
mented in  routine  NEODISP.  It  should  be  pointed  out  that  any 
modification  to  the  capabilities  of  the  translator  (i.e. 
being   able     to    display  other      types    of   quantities)  can   be 

made    by    modifying  these   routines   only. 
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LISTINGS 

This  appendix  contains  the  various  files  that  constituts 
the  translator  module.  The  submodula  that  displays  the  orig- 
inal geometry  is  listed  on  the  first  four  pages.  It  is 
imbedded  in  file  FRANK. FOR.  The  subaodule  that  takes  care  of 
the  post-processing  facilities  is  imbedded  in  file  DISP. FOR 
and  is  listed  in  the  remaining  pages.  It  has  been  mentioned 
previously  that  EPSA  had  been  slightly  modified  to  accomo- 
date color  graphics  capabilities.  The  only  interaction  of 
EPSA  with  the  translator  occurs  via  subroutine  calls.  All 
the  "calls"  occur  in  COMPUTE  (for  post-processing)  and 
REPORT  (for  the  original  geometry) .  A  labelled  COMMON 
called  FRANK  has  been  created  and  is  defined  in  the  routine 
AAA  as  requirad  by  EPSA.  The  regussts  for  PATRAN-G  outputs 
are  echoed  in  the  EPSA  output  file,  all  the  modifications 
for  that  purpose  having  been  made  in  routine  READIN.  Ths 
user  must  be  aware  that  the  size  of  blank  COMMON  array  A  has 
been  increased  to  store  the  deflections  in  x,  y,  z  direc- 
tions  instead  of   only  the  z    direction   previously. 
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Subroutine  orelim 

d  i  men  s  i  on  A  (  1  ) 

C  ommon  i  a ( 1  ) 

equivalence  (ia(l)»a(l)) 

common/ssi  ?e/  ibq( 1 ) »nq  j ,  nel tot#nlbd,nl oad,nbrect , 

1  nbouad,  isheet,  norofSi  nsotsi  nstrots,  nvots,  nhots, 

2  nsstvo»  nn  j  ,  nntot»  lqdso,  Maud*  lbcalc 
common  /cDara/  nsteos#  nbeg,  nend.  nsheet 

common  /stab  /  i bt (  1  )  r j ss i ze,  j soar ,   i ve I o» j st re* j xmas » j i e 1 m, j bmat , 

1  jllod,jloCD»jpret>jlhis,isrrnfiforc#jx1oc»jnqi#.inni,jnqbeg» 

2  j  1  s  i  d  e 


10 
11 


12 


common/ *rank/nfntot> 1 1  u 

CQMMON/TITRE/NTITLEC80) 

CHARACTER**?  BUFF 
CHARACTER*5  TITLE 
CHARACTER**  TIM 
CHARACTER*''  VER 

this  routine  comoutes  new  number  of  nodal  points 
for  any  given  sheet 


KJ=JMQ8 
nf n  t  O t  = 
1  1  t  voe  = 
1  1 kc  =  l 
1  1 i v  =  0 
1  1  id  =  0 
1  lnl=0 
Iln2=0 
1 1n3  =  0 
1 ln4  =  0 
1 ln5=0 
wri  ted 
forma  t ( 
wri  t  e ( 1 
format ( 


EG-JNNI 

nntot-IA(JNNI) 

25 


IA(JNQBEG-1 )-2*(KJ-2) 


lu, 10)  1 1  type,  1  Iid,1liv,llkc,llnl,1ln2,1ln3,lln4,lln5 

i2,5i8) 

lu,ll)  (NTITLE(I), 1=1,80) 

80  Al) 


tatting    care    of    second    oacket 

1 1 tyoe=26 

1 1 kc=l 

1 1 nl=nf ntot 

1 1 n2=nel tot 

writefllu,lO)  11  tvoe,  11  id,  11  i v» 1  1  kc» 11 nl , 11 nZ, 11 nl, 1 1n«,  1  InS 

call  date(8UFF) 

write(Uu,12)  BUFF 

format (A, 3x, ' 17:12:09' ,5X,' 1.4') 

return 

end 


subroutine  shconn 

dimension  A  (  1  ) 

common  i  a ( 1  ) 

equivalence  (ia(l)»a(l)) 

common/ssi  z  e  /  ibg( 1 ) ,  n  q  j ,  n  e  1 tot  #  n  1  b  d ,  n  1 oad,norect » 

1  noquad*  isheet.  norots,  nsotS»  nstrots,  nvpts,  nhpts# 

2  nsstvoi  nn i  ,    nntot»  lgdscw  lioud#  lbcalc 
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common  /coara/  a5tef>?»  nbeg,  nendi  nsheet 

common  /stab  /  ibt(l),jssize>jsoar,jvelo#jstre,jxmas,jielm,jbmat» 

1  jllcd»jlodD»joret»ilhis»jstrn,iforc»jxloc»jnqi,jnni,jnaoeg, 

2  j 1 s  i  de 

common/frank/nfntot/  1  1  u 

this  routine  will  write  the  element  corner  nodes  of  each 
element  on  neutral  file 

1  1  type  =  2 

1 1 kc  =  2 

1 1 nod=4 

1 1 i v=« 

LLN1=0 

LLN2=0 

LI_N3  =  0 

LLNttrO 

LLN5=0 

LLC0NF=0 

LLPI0=0 

LLCEI=0 

THET1=0. 

THET2=0. 

THET3=0. 

INITIALIZATION  OONE 

no  rev  =  0 

1  1 e 1 =  j  nn  i - j  nai 


do  200  k  =  t ,  1  lei 
1 1 row=I  A ( jnai  *k~\ ) 
number  of  elts  in  each  row 
do  100  j  =  1  , 1 1  row 
le1=IA(jnqeeg+K-l)+j-t 


nlel 1= 
nlel2= 
nlel 3= 
nl el  4  = 


+  no  rev 

+ 1 tnorev 

*  1  trior ev  +  1  1  row  *  1 

♦  norev  +  1  I  row  *  1 


ready  to  disolay  oacket 


1 1 id=LEL 

write(11u,80)  11 tyoe, 1  1  id, 1  1 iv,  1  1 kc / LLN 1 , LLN2 , LLN3, LLN4, LLN5 

80  format (i2,8i8) 

writed lu*81)     1 1 nod,LLCONF,LLPID,lLCEI , THET 1 , THET2, THET3 

81  f ormat ( i 8, 3i 8, 3el6.9) 

write(llu,82)  n1ell,nlel2,nle)3,nlel4 

82  format ( 10i8) 
1 00       cont i  nue 

norev=norev+l 1 row*  1 
200   continue 
return 
end 


subroutine  sheetf 
d  i  mens  ion  A ( 1 ) 
Common  i  a ( 1  ) 
eauivalence  (ia(l)»a(l)) 
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common/ssi ze/  ibq(l ) / no  j »nel tot  tnlblinloaa»nhr?cti 

1  noauaa»  isheet ,  norots*  nsots.  nstrctSi  nvots»  nnots* 

2  nsstvo/  n  n  j  ,  nntot,   I  }  d  s  o  >   l  i  q  u  d  »  l  b  c  a  l  c 
common  /CQara/  nsteosi  nbeq*  nend,  nsheet 

common  /stab  /  ibt(l)»jssize,jsoar,jvelo»jstre/jxmas,jielmfjbmat» 

1  Jl!oc1»j1odo>joret»jthis»jstrn,jforc>jx1oc>jnaif.jnni,jnabeq* 

2  j  1  s  i  "ie 

eommon/frank/nfntot* 1 lu 
character*l  qtyoe 

renumbers  the  nodes,  chanae  coordinates*  disolav  oacket  1 

1  tyoe=l 

1  id  =  0 

1  i  v  =  0 

1  kc  =  2 

lnl=0 

ln2  =  0 

ln3  =  0 

lna  =  0 

ln5  =  0 
nc  oun  =  0 

ldof=6 

1 i  c  f  =  1 
cjtyoe=,G' 

1  conf  =  0 

lei  d  =  0 
LSPC1=0 
LSPC2=0 
LSPC3=0 
LSPCU=0 
LSPC5=0 
LSPC6=0 

INITIALIZATION  DONE 

k  j  =  j  nabed- j  nn  i 
k  row  =  k  j -2 
do  200  j=l,krow 
I ood  on  new 


nb .  o  f  i"ow 


ncoun  =  ncoun-t-IA(  jnni  +j-l  ) 

11 ot  =  IA( j  nni  +  j -1 )-2 

do  100  1=1  ,  1  lot 
1 1 id=l 1 id+1 

xlloe=A(ncoun*2+jxloc*2*1 ) 
ylloc=A(ncoun*2fjx1oc*2*l+n 
skio  first  node  of  row 

zlloc=0. 

if  ( A( j vel o-2) .ne.0. )  then 

rcouro=l .  /  A  ( jvelo-2) 

theta=x 1 1 oc/fcouro 

xlloc=rcourb*sin(theta) 

z  1 loc=rcoupb*cos(theta)-rcourb 

endi  f 


if  (a( j vel o-l ) .ne.0. )  then 
rcouro=l  ./a( j  vel o- 1 ) 
the t  a  =  y  1  loc/rcouro 
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yl1oc=rcouro*sin(theta) 

z  1  loc=rcouro*cos(theta)-rcourb 

end i  f 

if((a(jvelo-l).ne.O.).and.(a(jvelo-2).ne.O.))  then 
wri te( 1  lu#  70) 
70  f ormat (' error , t wo  curvatures  are  non  zero') 

endi  f 

:       ready  to  disolay  oacket 

write(llu,30)  I  1 tyoe,  11  id#l  livr  1  lkc,  1  In  I,  1 1 n2# 1 1 n3» 1 1 n«r  IVnS 

80  format  (  i2,8i8) 

wri te( 1 1 u<81 )  t 1 1 oc / y 1 1 oc » 2 11 oc 

81  format (3el6.9) 

write(Hu»82)  1 H ef #gtyDe» 1 1 dof » 1 1 conf , 1  lei d,LSPCl rLSPC2# 
1   LSPC3  L5PC«,LSPC5,LSPC6 

82  format ( II ,al ,3i 8,2x,6i 1 ) 
1 00       cont  i  nue 

200    continue 
return 
end 


subroutine  enaneu 
COffimon/frann/nfntot# ' lu 
1  1  type=99 
1 1 id=0 
1 1 i  v  =  0 
1 1 kc=l 

write(llu»80)  11tyoe,llid»lliv,llkc 
80    format ( i  2, 3i  8) 
return 
end 
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N.ENO  ,     NSHEET  ,  N2BD 
NJOIN,  NRELAX,  ALPHA 


N3BD1 


SUBROUTINE  NEUDISP(DEFL,VELO) 

DIMENSION  A(l) 

COMMON  I  A ( 1  J 

EQUIVALENCE  (IA(1),  A(l)  ) 

DIMENSION  DEFL( 1 ) ,VELO( 1) 

COMMON  /CPARA/  NSTEPS  ,  NBEG  , 
1  N38D2,  INTRVL,  DELT,  NHTOT, 
5  LEN  S8X 

COMMON  /  SSIZE  /  I9GU),NQJ,NELT0T,N1BD,NL0AD,NBRECT, 

1  NBQUAD,  ISHEET,  NPRPTS,  NSPTS,  NSTRPTS,  NVPTS,  NHPTS, 

2  NSSTrP,  NNJ,  NNTOT,  LGDSP,  LIQUD,  LBCALC 

COMMON  /STAB  /  I8T( 1 ),JSSIZE, JSPAP, JVELO, JSTRE, JXMAS,JIELM, JB^AT, 

1  JL1BO, JLODP, JPRET,JLHI3, JSTRN, JFORC, JXLOC, JNQI , JNNI , JNQBEG , 
$    JLSIDE, JIELMCL, JSTIF, JDEFL, JFORLG, 

2  JIFPAR, JFLPAR, J  X COORD, JYCOORD, JDELTAX, JDELTAY, JVM A, JSEFX, 

3  JFLUFR,  JPRINC,  JVELRAD,  JGENFR , JPRES , JCSEP 

COMMON /CIO/  NIN,NOUT,NTHIST,NCORT,MCORT,NTPLOT,NVMA, ITITLEC20) 

COMMON/FRANK/NFNTOT,LLU,LLN,LLS,NDPLTS,NSPLTS 

C0MMON/TITRE/NTITLE(80) 

COMMON/CDELT/ISTEP.TIME 


3E 

:port 

BLANC 

2 

BLANC 

5 

BLANC 

a 

CPARA 

2 

CPARA 

3 

MSPARA 

3 

SSIZE 

2 

SSIZE 

3 

SSIZE 

a 

REPORT 

8 

STAB 

2 

STAB 

3 

STAB 

a 

STAB 

5 

STAB 

6 

10 


INTEGER    NSUBT1 (  80 ) , NSUBT2 ( 80 ) 

KJ=JNG8EG-JNNI 

NFNTOT  =  NNTOT-IA(JNNn-IA(JNQBEG-l  )-2*(KJ-2) 

LLU=19 

LLN=18 

LLM=17 

DEFMAX=0. 

MAXNO0=NFNTOT 

NWIDTH=5 

NDMAX=0 

LLID=0 

NCOUN=0 

KR0W=KJ-2 

DO  10  J=l,80 

NSU8T1 (J)=0 

NSUBT2(J)=0 


C  PUT  THE  TIME  STEP  ON  FIRST  SU8TITLE: 
C 

CLOSE(UNIT=LLM) 

OPE  N  (UNIT  =LLM,STATUS=' NEW) 

rtRITECLL^,*)  TI^E 

CLOSE  (UNIT=LLM) 

READ(LLM,99)  (NSUB T  1  ( J ) , J  =  1 , 80  ) 

CLOSE(UNIT=LLM) 

OPEN ( UNI T  =  LLM, STATUS =' OLD'  ) 

OPEN ( UN  I T=LLN,FORM=' UNFORMATTED' , ST A TUS=  '  NEW ' ) 


KROW  IS  THE  NEW  NUMB.  OF  ROWS  ,KJ  IS  THE  OLD 

NCOUM  COUNTS  NODAL  PTS  IN  EPSA,  LLID  COUNTS  FOR  ACTUAL  MODEL 

NUMR  COUNTS  NODAL  PTS  IN  EPSA 

FIRST  LOOP  TO  GET  MAX.  DEFORMATION  AND  NODE  NUMBER 

DO  200  J=t,KROW 
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NCOUNsNCOUN+IA(JNNI+J«l) 

LLPT=lA(JNNI+J)-2 

LLPT  IS  NUMB.  OF  POINTS  (FOR  PATRAN)  IN  EACH  ROW 

aE    TAKE  THE  RIGID  BODY  MOTION  OUT  BY  SUBTRACTING  THE  V  AND  ft 

DISPLACEMENTS  AT  END8L0CK  AT  EACH  NODAL  POINT 

NCOUN+2  IS  THE  POINT  N8 .  (EPSA)  FOR  ENDBLOCK 

RBYrOEFL(3*(NCOUNt-2)-l  ) 

RBZ=DEFL(3*(NCOUN+2) ) 

IF  (LIQUD.EQ.O)   THEN 

RBY=0. 

RBZ=0. 

ENOIF 

DO  100  L=ULLPT 

LLID=LLID*1 

NUMR=NCOUN+L+l 

DX=OEFL(3*NUMR-2) 

DY=DEFL(3*NUMR-1 )-PBY 

DZ=OEFL(3*NUM»)-R8Z 

IS  SHELL  CURVED, CHANGE  COORDINATES 

IF(  (ACJVELO-1 ) .NE.O.).OR. (A(JVELO-2) .NE.O.) )  THEN 

CALL  CHCOOPD(NUMR,DX,DY,DZ) 

ENDIF 
DD=AMAX1 CABS (OX), ABS(DY),ABS(DZ)) 
IF(ABS(DEFMAX)  .LT.DD)  THEN 
NOMAX=LLID 
DEFMflx=DD 

IF(ABS(AMIN1 (DX,DY,DZ) ) .EQ.DD)  THEN 
DEF^AX=-DEFMAX 
ENDIF 


:  WE  CHECKED   IF  DEFMAX  *AS  NEG. 

ENDIF 
100  CONTINUE 

200        CONTINUE 
I  OK  FOR  FIRST  CARD 

WRITE(LLN)  (NTITLEf  I) , I = 1 , 30 ) , NFNTOT , MAXNOD , DEFMAX , NDMAX , 
1  NWIDTH 
write(LLM,90)  TITLE,  NFM  TOT,  MAXNOD,  DEFMAX,  NOMAX,NWIDTH 

90  F0RMAT(A,2I8,E16.9,2I8) 

95    FORMAT(20A1, 218, £16.9,218) 

rtRITE(LLN)  (NSU8T 1 ( I ) , I = 1 , 30 ) 
WRITE(LLN)  (NSU8T2(I), 1=1,801 
WRITE(LLM,91 ) 

91  FORMAT(  'PINE  '  ) 
99    F0RMAT(80A1) 

LLID=0 
NCOUN=0 

I  SECOND  LOOP  TO  PICK  UP  DEFLECTION  AT  EACH  NODE  (OF  PATRAN  MODEL) 


DO  aOO  J=l,KROW 

NCOUN=NCOUN+IA(JNNI+ J-l ) 

LLPT=IA(JNNI+J)-2 

TAKE  OUT  RIGID  BODY  MOTION,  R8X,RBZ  DISPLACEMENTS  AT  END  BLOCKS 
ASIJMFD  TO  REPRESENT  RIGID  BODY  MOTIONS 
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RIGID  BODY  CONTRIBUTION 


NE.0.3)  THEN 


P8Y=0EFL(3*(NC0UN»2)-1 ) 
RBZ=OEFL(3*(NCO(JN  +  2)  ) 
IF  CLIQUD.EQ.O)  THEN 
RBY=0. 
RBZ=0. 
ENDIF 
C  IF  NO  FLUID  OPTION  DO  NOT  SU8TACT 
DO  300  L=1,LLPT 
LLID=LLID*1 
NUMR=NC0UN+L+1 
DX=OEFL(3*NUMR-2) 
DY=DEFL(3*NUMR-l)-RBY 
DZ=DEFL(3*NUMR)-R8Z 
DZLOC=DZ 
VELZ=VELO(3*NUMR) 

IFC(A(JVELO-l ) .NE.O.) .0R.(A(JVEL0-2) 
CALL  CHCOORO(NUMR,DX,DY,DZ) 
ENDIF 
WRITE (ULN)  LLID,DX,DY,DZ,DZLOC,VELZ 
*RITE(LLM,92)  LLID , DX , DY , OZ , DZLOC , VELZ 
300  CONTINUE 

aoo       CONTINUE 
C  CLOSE  FILE  OPENED  ON  UNIT  LLN 
CALL  CLOSE(LLN) 
92    F0RMAT(I8,5Elb.9) 
RETURN 
END 
C 

c 
c 

SUBROUTINE  CHCOORD ( N, X , Y , Z ) 

DIMENSION  A(l) 

COMMON  I  A  C 1 ) 

EQUIVALENCE  (IA(1),  A(l)  ) 

COMMON  /CPARA/  NSTEPS  ,  N3EG  ,  NEND  ,  NSHEET  ,  N2BD  ,  N3BD1  , 
1    N3BD2,  INTRVL,  OELT,  NHTOT,  NJOIN,  NRELAX,  ALPHA 
$  LEN  SBX 

COMMON  /  SSIZE  /  IBGU ),NQJ,NELT0T,N1BD,NL0AD,NBRECT, 

1  NBQUAD,  ISHEET,  NPRPTS,  NSPTS,  NSTRPTS,  NVPTS,  NHPTS, 

2  NSSTYP,  NNJ,  NNTOT,  LGDSP,  LIQ'JD,  L3CALC 
C 

COMMON  /STA8  /  IBT (l) » JSSIZE# JSPAR, JVELO, JSTRE, JXMASrJIELM, J8MAT, 

1  JL1BD, JL0DP,JPRET,JLHIS, JSTRN, JFORC, JXLOC , JNQ I , JNNI , JNQBEG, 
$    JLSIDE,JIELMCL, JSTIF, JDEFL , JFORLG, 

2  JIFPAR,JFLPAR,JXCOOPD,jYCOORD,JDELTAX,JDELTAY,JVMA,JSEFX, 

3  JFLUFR,  JPRINC,  JVELRAD,  JGENFR , JPRES , JCSEP 
COMMON/CIO/  NIN,NOUT,NTHI5T,NCORT,MCOPT,NTPL0T,NVMA,ITITLE(20) 
C0MWON/FRANK/NFNTOT,LLU,LLN,LLS,NDPLTS,NSPLTS 

DIMENSION  XMAT(3,3) 
DO  10  1=1,3 
DO  10  J=t ,3 
XMAT( I, J)=0. 
10   CONTINUE 
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IF  (ACJvELO-l).NE.O.J  THEN 
PC3URY=t ./A(JVELO-l  ) 
THETAl=A(JXL0C*2*N-l)/RC0URY 
XMATC |,l)sl, 
XMAT(2,2)=C0S(THETA1  ) 
XMAT(2,3)=SIN(THETA1  ) 
XMAT(3,3)=XMAT(2,2) 
XMAT(3,2)=-1 .*XMAT(2,3) 
CALL  PROD(XMAT,X, Y,Z) 
ENDIF 

IF  CURVATURE  IN  X  DIRECTION  IS  NON  ZERO: 


IF  (A(JVELO-2) .NE.O.)  THEM 
DO  20  1=1, 3 
DO  20  J=l,3 
XMATCI, J)=0. 
20     CONTINUE 

RC0URX=l./A(JVEL0-2) 

THETA2=A(JXL0Ct2*N-2)/RC0URX 

XMAT(2,2)=1. 

XMATd  ,  1  )=C0S(THETA2) 

XMAT(3,3)=XMAT(1, I) 

XMAT(3, 1)=-1 .*SIN(THETA2) 

XMATd,  3)=SIN(THETA2) 

CALL  PROD(XMAT,X, Y,Z) 

ENDIF 

RETURN 

END 


SUBROUTINE  PROD ( XWAT , X , Y , Z) 
DIMENSION  XMAT(3,3) 

THIS  ROUTINE  DOES  THE  MATRIX  PRODUCT  XMAT*(X,Y,Z) 

X1=X 

Y1=Y 

Z1  =  Z 

X=XMAT(l,l)*XltXMAT(l,2)*Yl+XMAT(l,3)*Zl 

Y  =  XMAT(2,1)*XH-XMAT(2,2)*Y1+XMAT(2,3)*Z1 

Z=XMAT(3, 1)*X1*XMAT(3,2)*Y1+XMAT(3,3)*Z1 

RETURN 

END 


SU8ROUTINE  PLOTDISPCDEFL, VELO) 

DIMENSION  DEFLC1 ),VELO( 1 ) 

COMMON/CDELT/ISTEP,TIME 

COMMON  /CPARA/  NSTEPS  ,  NBEG  ,  MEND  ,  NSHEET  ,  N2BD  ,  N3BD1  , 
1    N3BD2,  INTRVL,  DELT,  NHTOT,  NJOIN,  NRELAX,  ALPHA 
$  LEN  SBX 

COMMON/FRANK/NFNTOT,LLU,LLN,LLS,NDPLTS,NSPLTS 

CHECKS  IF  OUTPUT  FOR  PATRAN  IS  REQUESTED  BY  USER  IF  YES  CALL 
NEUDISP 

IF  (TIME. EQ. DELT)  KDISP=1 
TIMEI=FLOAT (NSTEPS) /FLOAT (NDPLTS) 
TI"EC=KDISP*TIMEI 
ITIME=JNINT(TIVEC) 


CPARA  2 
CPARA  3 
MSPARA  3 
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TIM£C=TIME/(ITI^E*OELT) 

ENDTIMsNSTEPS*DELT 

TESTaABS(TIM€C-JNINT(TIMEC) ) 

IF { (TEST. LE. 0.00001). OR. ( TIME. EQ.ENDTIM))  THEN 

CALL  NEUDISP(DEFL,VELO) 

KDISP=K0ISP+1 

ENDIF 

RETURN 

END 


SUBROUTINE  PLOTSTRE(STRE) 

DIMENSION  STRE(l) 

COMMON/CDELT/I STEP, TIME 

COMMON  /CPARA/  NSTEPS  ,  N8EG  ,  NENO  ,  NSHEET  ,  N2BD  ,  N3BD1 
1    N3BD2,  INTRVL,  DELT,  NHTOT,  NJOIN,  NRELAX,  ALPHA 
$  LEN  SBX 

COMMON /FRANK /NF NT OT , LLU , LLN , LLS , NDPL T S , NSPLTS 


CPARA  2 
CPARA  3 
MSPARA  3 


CHECKS  IF  OUTPUT  FOR  PATRAN  IS  REQUESTED  BY  USER  IF  YES  CALL 
NEUDISP 

IF  (TIME. EO. DELT)  KSTPE=1 

TIMEI=FLOAT(NSTEPS)/FLOAT(NSPLTS) 

TIMEC=KSTRE*TIMEI 

ITIME=JNINT(TIMEC) 

TIMEC=TIME/(ITIME*DELT) 

ENDTIM=NSTEPS*DELT 

TEST=ABS(TIMEC-JNINT(TIMEC) ) 

IF(  (TEST. LE. 0.00001). OR. (TI ME. EO.ENDTIM))  THEN 

CALL  NE'JSTRE(STRE) 

KSTRE=KSTRE+1 

ENDIF 

RETURN 

END 


SUBROUTINE  NEUSTRE ( STRE ) 

DIMENSION  A(l) 

COMMON  IA(1) 

EQUIVALENCE  (IAU),  A(l)  ) 

DIMENSION  DISP(10),STRA(10),DSTRE(10),DS(5),VEL(2) 

DIMENSION  STRE(l) 

COMMON  /CPARA/  NSTEPS  ,  NBEG  ,  NEND  ,  NSHEET  ,  N2BD  ,  N3BD1  , 
1    N38D2,  INTRVL,  DELT,  NHTOT,  NJOIN,  NRELAX,  ALPHA 
$  LEN  SBX 

COMMON  /  SSIZE  /  IBG(l) ,NQJ,NELT0T,N18D,NL0AD,NBRECT, 

1  NBQUAD,  ISHEET,  NPRPTS,  NSPTS,  NSTRPTS,  NVPTS,  NHPTS, 

2  NSSTYP,  NNJ,  NNTOT,  LGDSP,  LIQUD,  LBCALC 

COMMON  /STAB  /  IBT(1),JSSIZE, JSPAR, JVELO, JSTRE, JXMAS,JIELM,JBMAT, 

1  JL1BD, JLODP, JPRET, JLHIS, JSTRN, JFORC, JXLOC, JNQI, JNNI, JNQBEG, 
$    JLSIDE, JIELMCL/JSTIF, JDEFL , JFORLG, 

2  JIFPAR,JFLPAR,JXCOORD,JYCOORD,JDELTAX,JDELTAY,JVMA,JSEFX, 
J    JFLUFR,  JPRINC,  JVELRAD,  JGENFR , JPRE S , JCSEP 

COMMON/CIO/  NIN,NOUT,NTHIST,NCORT,MCORT,NTPLOT,NVMA,ITITLE(20) 
COMMON/FRANK /NFNTOT , LLU , LLN , LLS , MDPL T S , NSPLTS 
COMMON/TITRE/NTITLE(80) 
COMMON/CDELT/ISTEP, TIME 
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COMMON /SPAR  A /HEG( I ) ,E,GNU,.?HO,SGYLD,THICK,CURV(2) 

INTEGER   NSUBTt  (80)  ,NSU9T2(80) 

N*IDTH=3l 

LLS  =  t6 

LLZ=tS 

00  10  J=l , 30 
tO       NSU8T2(J)=0 
C  TAKE  CARE  OF  THE  2N0  TITLE  ,  INOICAT.  TIME  STEP 
C 

CLOSE(UNIT=LLZ) 

OPEN (UNITsLLZ, STATUS* • NEW ) 

wRITECLLZ,*)  TIME 

CLOSE(UNIT=LLZ) 

REA0(LLZ,801)  (NSUBT1 (I) ,1=1,80) 

CL0SE(UNIT=LLZ) 

OPEN(UNIT=LLZ,STATUS='OLD'  ) 
801   FORMATC80A1) 

OPEN ( UN  I T  =  LLS,FORM=' UNFORMATTED" , ST ATUS= ' NEW ' ) 
C 

c 

C  WE  ARE  GOING  TO  DISPLAY  THE  FIRST  THREE  RECOROS  (TITLES) 
C 

WRITE(LLS)  (NTITLE(I),T=1.80) ,NWIDTH 
WRITE(LLS)  (NSU8T1(I),I=1,80) 
WRITE(LLS)  (NSUBT2(I),I=1,80) 
C 

C         WE  ARE'  GOING  TO  SCAN  ALL  ELEMENTS  AS  IN  SHCONN, 
C         GET  THE  STRESSES  PERFORM  ROTATION  IF  NECESSARY 
C         LLEL:  NUM.  OF  ROWS  ELTS;  LLID:ELT  NUMBER 
C         NPREV:  NODE  NUM. (FOR  EPSA)  OF  LAST  POINT  OF  ROW  JUST 
C         BELOW  THE  ELEMENT  ROW  X;  LLROW:  NUM.  OF  ELTS  IN  EACH  ROW 
C 

NPREV=0 
C  IA  ELTS  *t(FOR  NODES)  +2  ARTIFICIAL  NOOES  FOR  NPREV 

NSHAPE=a 
LLEL=JNNI-JNQI 
00  200  K=1,LLEL 
C 

LLR0W=IA(JNQI+K-1) 
NPREV=NPREV+LLR0W*3 
C 

DO  100  J=t, LLROW 
LLID=IA(JNQBEG+K-1)*J-1 
DS(1)=STRE(9»(LLID-!)*1) 
DS(2)=STREC9*(LLID-l)t2) 
DS(3)=0. 
C  DSU,5)  WILL  BE  THE  STRESSES  IN  LOCAL  COORD. 
DS(U)=0S(1) 
DS(5)=0S(2) 
C  VON  MISES  AT  TOP  AND  BOTTOM  OF  SHELL: 
T0XY=STRE(9*(LLID-l)+3) 
DSX=STRE(9*(LLID-l)+«) 
DSY=STRE(9*(LLlD-t )+5) 
DSX=DSX*6. /THICK 
DSY=DSY*6. /THICK 
DSX1=0S(1)+DSX 
DSY1=0S(2)+0SY 
DSX2  =  DS(  D-OSX 
DSY2=OS(2)-OSY 
C 
C   VON  MISES  CALCULATIONS  MODIFIED  AFTERWARDS  FOR  MEMB .  1-BENDING  STRESS 
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C   NO  SHEAR  AT  TOP  OF  SHELL, GET  STRESSES  AT  TOP  ANO  BOTTOM  AMO  TAKE 
C   MAX.  VALUE 
C 

C  TOP: 

SIGMt=OSXl 
SIGM2=OSYl 

V0NM1=SQRTCSIGM1*SIGM1-SIGM1*SIGM2+SIGM2*SIGM2) 
C  BOTTOM: 

SIGMt=OSX2 
SIGM2=OSY2 

V0NM2=SQRT(SIGM1*SIGM1-SIGM1*SIGM2+SIGM2*SIGM2) 
C 

VONMsAMAXl (VONM1 , V0NM2) 
C     O.K.  FOR  VON  MISES  (VONM) 

C     IS  THE  SHELL  CURVED,  NLELI'S  ARE  NOOES  SURROUNOING  ELT.  (EPSA) 
C 

IF(CA(JVELO-1).NE.O.) .0R.(A(JVEL0-2) .NE.O.))  THEN 
NLEL1=J+NPREV+1 
NLEL2=J>NPREV+2 
NLEL3=J+ 1 +NPREV fLLROW* 1 +2 
NLELU=J>NPREV+LLR0W'H+2 

X=(A(JXLOC+2*NLELl-2)+A(JXLOC  +  2*NLEL2-2))/<4. 
X=X+(A(JXLOC+2*NLEL3-2)+A(JXLOC+2*NLEL«-2) )/U. 
Y=(A(JXLOC+2*NLELl-l)*A(JXLOC+2*NLEL2-n  )/«. 
Y=Y*(A(JXL0C+2*NLEL3-1  )  ♦ A ( JXL0C*2*NLEL4- I ) ) /4 . 
CALL  CHELEM(X, Y,DS) 
ENOIF 
C 

C  READY    TO    DISPLAY    RECORD 

C 

rfRITE(LLS)    LLID,NSHAPE, (D I SP ( I ) , 1= I , 1 0 ) , 

1  (STRA(I) ,1=1, 10),DSTRE(1),DS(5) ,DS(6), 

2  DSTREm,  (DS(I),I  =  1,3), 

3  (DSTRE(I) ,1=8,10) , VONM 

C  WRITEC15,80)  LLID.NSHAPE, (OS ( I ) , 1  =  1 , 5) , VONM 

100      CONTINUE 
200   CONTINUE 

CALL  CLOSE(LLS) 
80    F0RMAT(2ia,6E15.8) 

RETURN 

END 


SUBROUTINE  CHEL£M( X , Y , OS ) 

DIMENSION  A(l) 

COMMON  I A  (  1  ) 

EQUIVALENCE  CIA(l),  A(l)  ) 

COMMON  /CPARA/  NSTEPS  ,  NBEG  ,  NEND  ,  NSHEET  ,  N2BD  ,  N3BD1  , 
1    N38D2,  INTRVL,  DELT,  NHTOT,  NJOIN,  NRELAX,  ALPHA 
S  LEN  SBX 

COMMON  /  SSI2E  /  IBG(l) ,NQJ,NELT0T,Nl8D,NL0A0,N6RECT, 

1  N8QUAD,  ISHEET,  NPRPTS,  NSPTS,  NSTRPTS,  NVPTS,  NHPTS» 

2  NSSTYP,  NNJ,  NNTOT,  LGDSP,  LIQUD,  LBCALC 

COMMON  /STAB  /  IBT(1),JSSIZE, JSPAR,JVELO, JSTRE, JXMAS , JIELM, JBMAT , 

1  JL18D,JL00P,JPRET,JLHIS,JSTRN,JF0RC,JXL0C,JNOI,JNNI,JNQBEG# 
S    JL5IDE,JIELMCL,JSTIF,JDEFL,JF0RLG, 

2  JIFPAR,JFLPAR,JXCOORD,JYCOORD,JDELTAX,JOELTAY,JVMA,JSEFX, 

3  JFLUFR,  JPRINC,  JVELRAD,  JGENFR , JPRES, JCSEP 
DIMENSION  XMAT(3,3) 

DIMENSION    DS(3) 
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00  10  1  =  1  ,  3 
00  10  Jsl,3 
XMAT{I,J}sO. 

10   CONTINUE 

THIS  ROUTINE  CHANGES  THE  STRESSE  IN  AN  ELEMENT  INTO  A 
IN  A  GLOBAL  RECTANGULAR  SYSTEM  .GIVEN  LOCAL  COORD.  OF 
CENTROID  X,r 

IF  CURVATURE  IN  Y  OIR.  IS  NON  ZERO  , CALCULATE 
THE  ROTATION  MATRIX  AT  EACH  POINT 
IF  (A(JVEL0-1 ) .NE.O. )  THEN 
RC0URY=1 ./A(JVEL0-1 ) 
THETA1=Y/RC0URY 
XMAT(1, 1 )=1 . 
XMAT(2,2)=C0S(THETA1 ) 
XMAT(2,3)=SIN(THETA1 ) 
XMAT(3,3)=XMAT(2,2) 
XMAT(3,2)=-1 .*XMAT(2,3) 
CALL  PROD(XMAT,DSU ) , OS ( 2 ) , OS ( 3 ) ) 
ENOIF 

:  IF  CURVATURE  IN  X  DIRECTION  IS  NON  ZERO: 
IF  CA( JVELO-2) .NE.O.)  THEN 

DO  20  1=1 ,3 

DO  20  J=l,3 

XMAT(I,J)=0. 
20     CONTINUE 

RC0URX=1 ./A(JVEL0-2) 

THETA2=X/RC0URX 

XMAT(2,2)=1 . 

XMATC1, 1)=C0S(THETA2) 

XMAT(3,3)=XMAT(1 , 1) 

XMAT(3,1)=-1.*SIN(THETA2) 

XMAT(1,3)=SIN(THETA2) 

CALL  PROO(XyAT,OS(l) , OS ( 2 ) , DS ( 3) ) 
ENOIF 
RETURN 
END 
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